(navigation image)
Home American Libraries | Canadian Libraries | Universal Library | Community Texts | Project Gutenberg | Children's Library | Biodiversity Heritage Library | Additional Collections
Search: Advanced Search
Anonymous User (login or join us)
Upload
See other formats

Full text of "Naval research logistics quarterly"

1 1 ,< *>s 



NfiVflL RESEARCH 

LOGISTICS 
QUflRTERLy 



JUNE 1970 

VOL. 17, NO. 2 fe 

re 



13 OCT 7 0% 




OFFICE OF NAVAL RESEARCH 



NAVSO P-1278 



NAVAL RESEARCH LOGISTICS QUARTERLY 



EDITORS 

Rear Admiral H. E. Eccles, USN (Retired) 
The George Washington University 

F. D. Rigby 
Texas Technological College 

S. M. Selig 

Managing Editor 

Office of Naval Research 

Arlington, Va. 22217 



O. Morgenstern 
Princeton University 

D. M. Gilford 
U.S. Office of Education 



ASSOCIATE EDITORS 



R. Bellman, RAND Corporation 

J. C. Busby, Jr., Captain, SC, USN (Retired) 

W. W. Cooper, Carnegie Mellon University 

J. G. Dean, Captain, SC, USN 

G. Dyer, Vice Admiral, USN (Retired) 

P. L. Folsom, Captain, USN (Retired) 

M A. Geisler, RNAD Corporation 

A. J. Hoffman, International Business 

Machines Corporation 
H. P. Jones, Commander, SC, USN (Retired) 
S. Karlin, Stanford University 
H. W. Kuhn, Princeton University 
J. Laderman, Office of Naval Research 
R. J. Lundegard, Office of Naval Research 
W. H. Marlow, The George Washington University 
B.J. McDonald, Office of Naval Research 
R. E. McShane, Vice Admiral, USN (Retired) 
W. F. Millson, Captain, SC, USN 
H. D. Moore, Captain, SC, USN (Retired) 



M. i. Rosenberg, Captain, USN (Retired) 

D. Rosenblatt, National Bureau of Standards 

J. V. Rosapepe, Commander, SC, USN (Retired) 
T. L. Saaty, U.S. Arms Control and 
Disarmament Agency 

E. K. Scofield, Captain, SC, USN (Retired) 
M. W. Shelly, University of Kansas 

J. R. Simpson, Office of Naval Research 

J. S. Skoczylas, Colonel, USMC 

S. R. Smith, Naval Research Laboratory 

H. Solomon, The George Washington University 

I. Stakgold, Northwestern University 

E. D. Stanley, Jr., Rear Admiral, USN (Retired) 

C. Stein, Jr., Captain, SC, USN (Retired) 

R. M. Thrall, University of Michigan 

T. C. Varley, Office of Naval Research 

C. B. Tompkins, University of California 

J. F. Tynan, Commander, SC, USN (Retired) 

J. D. Wilkes, Department of Defense, OASD (ISA) 



The Naval Research Logistics Quarterly is devoted to the dissemination of scientific information in logistics and 
will publish research and expository papers, including those in certain areas of mathematics, statistics, and economics, 
relevant to the over-all effort to improve the efficiency and effectiveness of logistics operations. 

Information for Contributors is indicated on inside back cover. 






The Naval Research Logistics Quarterly is published by the Office of Naval Research in the months of March, June, 
September, and December and can be purchased from the Superintendent of Documents, U.S. Government Printing 
Office, Washington, D.C. 20402. Subscription Price: S550 a year in the U.S. and Canada, S7.00 elsewhere. Cost of 
individual issues may be obtained from the Superintendent of Documents. 

The views and opinions expressed in this quarterly are those ofthe authors and not necessarily those of the Office 

of Naval Research. 

Issuance of this periodical approved in accordance with Department of the Navy Publications and Printing Regulations, 

NAVEXOS P-35 



Permission has been granted to use the copyrighted material appearing in this publication. 



MARKOV CHAIN ANALYSIS OF A SITUATION WHERE CANNIBALIZA- 
TION IS THE ONLY REPAIR ACTIVITY 



Alan J. Rolfe 



The Rand Corporation 
Santa Monica, California 



ABSTRACT 



This paper considers a group of S identical aircraft, each of which is partitioned into 
K parts which fail exponentially. The only way in which a failed aircraft can be repaired is by 
cannibalizing its out-of-commission parts from other failed aircraft. The evolution of the 
number of good aircraft over time is governed by the transient behavior of an absorbing 
Markov chain. We can therefore study this behavior by matrix multiplication although the 
computational problem grows large for K 5= 3. Some numerical results and some approxi- 
mations are also provided. 

1. INTRODUCTION 

This paper presents a model of a system subject to cannibalization. The system is a simple special 
case of the very general coherent structure considered by Hirsch, Meisner, and Boll (HMB) in [3]. 
It differs, however, from the model of HMB in its treatment of time. 

In Sec. 2 the model is described in terms of a specific example. The reader can no doubt imagine a 
number of other situations where the same model would be appropriate. We also make clear how our 
model treats time differently from the model of HMB. Section 3 presents the Markov chain transition 
matrix which characterizes our problem, while Sec. 4 explores the computational chore associated 
with using the model. In Sec. 5 and 6, we exhibit several modest numerical examples and suggest 
means of lessening the computational chore for larger problems. Section 7 contains a special technique 
which is feasible for small-scale problems. 

Although this paper does not completely and explicity solve the problem it poses, we hope that 
the results given will provide quantitative insights into the effects of cannibalization which will lead 
to further work in the field. 

2. THE PROBLEM 

We consider a group of S identical aircraft operating under the following conditions: 

(a) Every aircraft includes K parts or assemblies. The time-to-failure of the ith part is exponentially 
distributed with distribution function, 1 — e~V. Futhermore, the K parts are assumed to be noncritical 
items so that no matter how many of the K assemblies fail during a mission it is always possible for an 
aircraft to return safely to base. 

(b) A failed part can neither be repaired nor resupplied. However, there is sufficient repair capa- 
bility to cannibalize unfailed parts from other failed aircraft in an attempt to minimize the number of 
out-of-commission aircraft. We assume that initial supplies of spare parts are zero, although this 
condition can be relaxed as we indicate later. 

(c) The aircraft fly missions of fixed length T in which every "good" aircraft flies every mission. 
We assume that there is enough time between consecutive mission takeoff times for all aircraft to 

151 



152 A. J. ROLFE 

return from the previous mission and for all cannibalizations to be performed. We indicate below how 
the "every good aircraft flies every mission" rule can be modified. 

The basic quantity which the analysis below will provide is the expected number of good aircraft 
after t missions, t=\, 2, . . ., denoted by E[A{t)]. This information will, for example, allow us to 
measure the changes in force capability over time which cannibalization affords under the conditions 
above. 

It is important to distinguish the way in which our model differs in its treatment of time from that 
of HMB. In our case, only those parts which happen to be in the air during a given period of time are 
subject to risk. Good parts lying in partially cannibalized aircraft on the ground are not subject to 
failure. In effect, each part carries its own clock which is turned on when the part is sent up into the 
air and turned off when the part is on the ground. In the paper of HMB, however, every part that has 
not previously failed is subject to risk. There is a single clock which is turned on at time 0, and while 
it runs all parts are subject to failure. In Sec. 6 below we use the above fact to identify one of the 
results of HMB as a lower bound for E[A(t)]. 

3. THE MARKOV TRANSITION MATRIX 

If we let rii(t) be the number of copies of part i which survive t missions, then n(t) = (ni(t), . . ., 
nA-(O) represents the number of parts of each type which survive t missions. This is the basic state 
vector of our process. As n.i{t) is an integer and ^ rii{t) ^ S for i=l, 2, . . ., K, the process has 
(S + 1) K states for each t. However, for our purposes we can lump all states having one or more ni{t) =0 
into a single "zero" state, thereby reducing the number of states to S K + 1. It is readily seen that the 
process is Markovian and that its evolution from mission t to t+l can be described by a one-step 
Markov chain transition matrix P of dimension S K + 1 by S K + 1. If the rows and columns are ordered 
lexicographically, P will be a lower triangular matrix with 0's above the principal diagonal. 

It is convenient to label row and column m by d m , where 

d\ =(zero state); 
<k =(1, . . ., 1); 



d S K +1 =(S, . . .,S). 
Then the matrix element py, where dt= (t'i, . . ., i/<), dj= (j\ . . .,Jk), and M,= min ii , is given by 



i=\ 



K / M. \ 

(2a) f[ ' ,„ ) (e- K i T ) j i- i ' +M i(l-e- K i T ) i i-Ji i 

jJi \ji — ii + Mj/ 

if i,y>l and;, =S i, <;, + M„ for all /=1, . . .,K; 

(2b) if i, j > 1 but ii < ji or i t > ji + Mi for at least one / = 1 , . . . , K; 

(2c) lifi=j = l; 

(2d) Oif £=1,j>1; 



MARKOV CHAIN CANNIBALIZATION ANALYSIS 153 

i 
(2e) 1~2 Pidi£i>l,j=l- 



rf=2 



These entries become apparent if we realize that Mi copies of each part are tested on each mission 
and that each part is assumed to fail independently of other parts. Once the zero state is reached the 
"experiment" ends as all aircraft are out-of-commission. 

4. THE COMPUTATIONAL CHORE 

If we let ri(<) be the row vector of unconditional state probabilities after t missions, then from 
Markov chain theory we know that 

(3) n(*)=n(0)F, «=i,2, . . ., 

where 11(0) = (0, . . .,0,1) as the process starts in state (S, S. . . . , S) and the states of II() are 
ordered lexicographically as are those of P. Furthermore, P represents an absorbing Markov chain 
[4] with d\ the absorbing state and all other states transient. From (3) we see that to go from Tl(t) to 
Yl(t-\-\) we must multiply 11(0 , a 1 byS K + 1 matrix, byP, anS^-f 1 byS A + 1 matrix. This operation re- 
quires roughly (S 2K /2) multiplications and the same number of additions. The major problem in using 
the model is one of computational feasibility as (S 2A /2) can be a large number. For example, if S = 24 
and K = 3, (S 2K /2)= 1.5 X 10 8 . Thus it appears that either S or K must be fairly severely restricted 
if we are to perform the computation in an acceptable amount of time; the following sections describe 
some numerical results which suggest ways of lessening the computational chore. 

The first two moments of the number of good aircraft after t missions are given by 



(4) E[A(t)]= S 2 nj(t)Mj 
and j=1 

(5) E[A(t)*] = S! 2 Uj(t)Ml 

where 

Mj= min {ji\dj={ji, . . ., >*)}• 

1=1, . . .,K 

Since Mj is independent of t for each;, the computation of E[A{t)] should proceed rapidly once 
Yl(t) has been determined. As the sole absorbing state d\ of the Markov process has Mi = it is clear 
from (4) that lim E[A(t)]=0; this is of course what one would expect given the no-repair situation 
being modeled. The computation of 11(f) would generally be iterated until E[A(t) ] was fairly close to 
zero. 

5. NUMERICAL RESULTS FOR K = 2 

In order to simplify the computational task, we first limit ourselves to the case of two cannibalizable 
parts per aircraft. There remain, then, two dimensions to the problem: S, the number of aircraft under 
test, and g = {e~ x,T , e~ k2T ), the vector whose components are the probabilities of nonfailure of each 
part on one mission of length T. Table 1 presents results for q fixed and S variable, while in Table 2 
q is variable and S fixed. In both cases the quantity of interest is E[A(t)], the expected number of 
good aircraft after t missions. 



154 



A. J. ROLFE 



Table 1 



\ s 

t \. 


Any S, no 
cannibalization; 

(0.72)' 


2 


5 


10 


13 


1 


0.720 


0.734 


0.762 


0.777 


0.783 


2 


0.518 


0.557 


0.600 


0.622 


0.628 


3 


0.373 


0.431 


0.480 


0.500 


0.504 


4 


0.268 


0.336 


0.384 


0.402 


0.405 


5 


0.193 


0.263 


0.308 


0.323 


0.326 


6 


0.139 


0.206 


0.247 


0.260 


0.261 


7 


0.100 


0.161 


0.198 


0.209 


0.210 


8 


0.072 


0.125 


0.159 


0.168 


0.169 


9 


0.054 


0.097 


0.128 


0.135 


0.136 


10 


0.037 


0.075 


0.102 


0.108 


0.109 


11 


0.027 


0.058 


0.082 


0.087 


0.087 


12 


0.019 


0.044 


0.066 


0.070 


0.070 


13 


0.014 


0.034 


0.053 


0.056 


0.056 


14 


0.010 


0.026 


0.042 


0.045 


0.045 


15 


0.007 


0.019 


0.033 


0.036 


0.036 



(Entries are E[A(t)]IS; q = (0.9, 0.8). For the no cannibalization case 
E[A(t)]/S = (0.72)' results from the fact that 0.72 is the probability that 
one aircraft survives one mission and that we have a sequence of t bi- 
nomial trials.) 



From Table 1 we see that for a fixed value of t, the ability to cannibalize produces an increase in 
E[A(t)]IS which also increases with S. However, we note that as S increases the quantity E[A(t)]/S 



Table 2 



t \v 


0.9, 0.8 


0.8, 0.8 


0.7, 0.7 


0.6, 0.6 


0.8, 0.5 


0.5, 0.5 


0.9, 0.4 


0.5, 0.4 


1 

2 

3 

4 

5 

6 

8 

10 

12 

14 


0.777 
0.804 
0.803 
0.806 
0.805 
0.802 
0.803 
0.804 
0.804 
0.804 


0.731 
0.768 
0.773 
0.775 
0.777 
0.776 
0.776 
0.775 
0.773 


0.620 
0.661 
0.667 
0.668 
0.668 
0.667 
0.660 


0.514 
0.558 
0.563 
0.558 
0.556 
0.550 


0.494 
0.506 
0.504 
0.503 
0.503 


0.412 
0.456 
0.456 
0.450 


0.399 
0.400 
0.401 
0.400 


0.353 
0.391 
0.389 
0.382 

































































(Entries are E[A(t)]/E[A(t- 1)]; S = 10.) 



converges for values of S around 10. The entries for S — 13 are not much different from those forS = 10. 
Hence it seems that we can learn as much about E[A(t)]/S by studying 13 aircraft as by studying 24 or 
48 and the computational job is markedly easier. Of course, this conclusion is strictly valid only for 
K — 2 and q = (0.9, 0.8). A question of interest (as yet unanswered) is, "At which value of S does con- 
vergence occur, given a particular K and q?" 

In Table 2 we exhibit the behavior of successive values of E[A(t)] for a group of 10 aircraft and a 
variety of g-vectors. One fact is immediately apparent: For t 2= 2, E[A(t)]/E[A(t— 1)] =p, p a constant. 



MARKOV CHAIN CANNIBALIZATION ANALYSIS 155 

This is equivalent to saying £'[/i(f)]/S = p'" 1 , t Ss 2, which implies 

(6) E[A(t)]IS = p ( P y-\t^l, 

where 

pomEtA(l)]IS. 

This provides an immediate analogy to the no-cannibalization case, where we know that E[A{t)]/S = A.', 
«Sl, with A = <7i<7 2 where q = (qi , q 2 ) . (It was found that for S = 5 the geometric behavior was less stable 
than that characterized by Table 2 for S= 10). 

The next section will suggest several possible ways around the computational problem. 

6. APPROXIMATIONS FOR K>2 

In Sec. 5 we indicated that for K= 2, any S 3= 10 and any q we could approximate the behavior of 
E[A(t)]/S by studying that of the case S— 10. A typical computer run on the 360/65 for S= 10 and 
q= (0.9, 0.8) took approximately 7 sec. If K were increased to 3, the information in Sec. 4 would 
suggest a run time of roughly 10 min., still a feasible amount of time. Hence for K= 2,3 we can study 
a group of 10 aircraft and be fairly confident that the E[A (t)]IS so obtained closely approximates that 
for any larger value of S and the same q. 

For K 25 4, however, studying even 10 aircraft becomes prohibitively costly as each time K is 
increased by 1, computation times are increased by a factor of 100. Thus, although the geometric 
behavior exhibited by E[A(t)]/S in (6) would probably still hold, the problem of determining p and p 
could no longer be handled by numerical means. 

An approximation for E[A(t)] is given by (7), with q= {q u ■ ■ • , <7a ) and a=mjn </,, 

(7, 5f*^,„i 

Relation (7) is based on the imprecise reasoning that A(t) is equal to the number of survivors of the 
part-type for which qi is minimized. 

The numbers in Table 2 suggest that the RHS of (7) is often an upper bound for E[A(t)]IS rather 
than an approximation; however, column 5 of that table shows that this is not always so. For t=l, the 
relationship E[A( 1)]/S =S a follows readily from min Xj^Xi, i=l, . . . , K, 

i = 1 K 

K 

A lower bound for E[A(t)]/S is given by (8), with |8 = r[ qi, 

The referee suggested the following better lower bound, which can be derived from formulae 4.3 
and 4.9 of [3]. This result can also be derived directly without reference to [3] and is based on the 
supposition that every surviving part be tested on every mission. 



(9) E[A(t)] 






The referee provided the results shown in Table 3 whose entries are the lower bound, (9), for the case 
shown in Table 1. Expression (9) appears to be an excellent lower bound in this instance. 



156 



A. J. ROLFE 
Table 3 



\ s 










t \^ 


2 


5 


10 


13 


1 


0.734 


0.759 


0.777 


0.782 


2 


0.554 


0.598 


0.620 


0.626 


3 


0.423 


0.473 ■ 


0.495 


0.501 


4 


0.323 


0.374 


0.396 


0.401 


5 


0.247 


0.296 


0.316 


0.320 


6 


0.187 


0.233 


0.252 


0.256 


7 


0.142 


0.184 


0.201 


0.204 


8 


0.106 


0.144 


0.160 


0.163 


9 


0.080 


0.113 


0.127 


0.130 


10 


0.059 


0.088 


0.101 


0.103 


11 


0.044 


0.068 


0.080 


0.082 


12 


0.032 


0.052 


0.063 


0.065 


13 


0.024 


0.040 


0.049 


0.052 


14 


0.017 


0.030 


0.039 


0.041 


15 


0.013 


0.023 


0.030 


0.032 



7. A RECURSIVE RELATION FOR THE GENERATING FUNCTION OF THE STATE 
PRORABILITIES 

The fact that the P-matrix, as given by expression (2), is lower triangular allows us to develop a 
recursive relation which, at least in theory, produces a closed form expression for the generating func- 
tion of the state probability vector, Yl{t). Inverting these generating functions presents a barrier if 
S and Kare large; however, we present the method as it might sometimes be desirable to study a problem 
with small S and K. 

Setting N=S K +l, we let {ptj, i, j=l, . . ., N} be the elements of the P-matrix. If II(t) = (IIi(f), 
. . ., Il/v(t)) is the vector of state probabilities, we define the generating functions 



(10) Ui(Z) = £ Yli(t)Z', i = l, 2, . . ., N. 

t = 

The basic Markov chain relation is 



(11) 



nu) = n{t-\)P. 



Using (10), (11), and the fact that pjj = 0,j> i, we can derive 

U N (Z) =p NN ZU N (Z) + 1 



This is equivalent to 



(12a) 
(12b) 



Yl i (Z) = ^ Pji Zn j (Z), i=N-l,N-2, . . ., 1. 



n N {Z) = \j{\-p NN z) 



n,(Z) = (z 2 p ji n j (Z))i(i-p ji z), i=N-i,N-2, . . ., 1. 

j = i + l 



MARKOV CHAIN CANNIBALIZATION ANALYSIS 



157 



From (12a) we know at once that 



Us(t) = (p NN y, t^O. 



Then (12b) can be solved recursively as IT(Z) depends only upon flj(Z), j>i. The difficulty in a 
large problem is that inverting Yli(Z) is not trivial as it was in (12a) because the generating function will 
have a large number of multiple roots. There do exist both exact and approximate methods for such 
problems, some of which are found in [2, Chap. XI]. We present a short numerical example to illustrate 
the technique: 

S =2 

K = 2 

q= (0.5, 0.4) 



From expression (2), we find that 



P = 



1 














0.8 


0.2 











0.5 


0.3 


0.2 








0.6 


0.2 





0.2 





_0.52 


0.24 


0.08 


0.12 


0.04 



Then expression (12) produces 

n 5 (Z) = i/u-o.04Z) 

fl 4 (Z) = (0.12Z)/(1 -0.2Z) (1 -0.04Z) 

n 3 (Z) = (0.8Z)/(1 -0.2Z) (1 -0.04Z) 

n 2 (Z) = (0.24Z)/(l-0.2Z) 2 (l-0.04Z) 

n I (Z) = (0.0016ZM5 + Z)(65-Z)/(l-Z)(l-0.2Z) 2 U-0.04Z). 

If we use methods in [1] to invert these generating functions and then use the relation 



we find 



E[AU)] = o{n l u)} + i{n 2 (t) + iu(t) + n 4 (t)} + 2{Ti i (t)}. 



^(*)]=|«(.2V+|(.2)'+|(.04)',«s»0. 



Thus we have found a closed form expression for E[A(t)]. 



8. VARIATIONS OF THE MODEL 



Several of the assumptions made in Sec. 2 can be relaxed without materially affecting the develop- 
ment outlined above. The only change is in the Markov transition matrix P: 



158 A. J. ROLFE 

(a) There may be an initial supply, s,, of spares for part i, i= 1, 2, . . . , K. In this case the possible 
states of part i are 0, 1, . . .,5, + S and the P-matrix has 

f[( Si + S) + l 

rows and columns. This would make the computational task mentioned in Sec. 4 more difficult; 

(b) The assumption that "every good aircraft flies" every mission" can be replaced by the re- 
quirement that only Q aircraft, if available, fly every mission. This would be accomplished by sub- 
stituting Q for Mi in expression (2a) whenever Q < Mi. If this were done, the P-matrix would be made 
more sparse while E[A(t)] would approach zero less rapidly because the Markov process would take 
more time to work its way from the initial state to the absorbing state. 

9. ACKNOWLEDGMENTS 

The author would like to thank Dr. D. M. Landi of The Rand Corporation for his part in bringing 
the model to the author's attention. The author is especially indebted to an anonymous referee who 
supplied many useful comments, all of which are incorporated in this paper. 



REFERENCES 

[1] D'Azzo,J. J. and C. H. Houpis, Feedback Control System Analysis and Synthesis (McGraw-Hill, New York, 1960). 

[2] Feller, W., An Introduction to Probability Theory and Its Applications (Wiley, New York, 1957), Vol. 1. 

[3] Hirsch, W., M. Meisner, and C. Boll, "Cannibalization in Multicomponent Systems and the Theory of Reliability," Na'v. 

Res. Log. Quart. 7.5, 335-359 (1968). 
[4] Kemeny, J. G. and J. L. Snell, Finite Markov Chains (Van Nostrand, Princeton, New Jersey, 1960). 



SOME ESTIMATES OF RELIABILITY USING INTERFERENCE 

THEORY 



M. Mazumdar 

Westinghouse Ktscarch Laboratories . 
Pittsburgh, Pennsylvania 15235 



ABSTRACT 

According to "interference theory" of reliability, a component fails if the maximum 
stress exceeds the component's strength. Assuming that both these quantities are random 
and their distributions are normal, we obtain in this paper some point and interval estimates 
of reliability when the stress distribution is known and a few observations exist on com- 
ponent strengths. 

1. INTRODUCTION 

In recent years, many writers have attempted to analyze the reliability of a component by the 
application of probability arguments to a certain physical model of failure. According to this model 
(which has come to be recently known in the literature as "interference theory"), a component fails 
if at any moment the applied stress (or load) is greater than its strength or resistance. The stress is a 
function of the environment to which the component is subjected, and its value at any point of time 
is considered to be a random variable. The strength of the component is measured by the stress re- 
quired to cause failure. Strength depends on material properties, manufacturing procedures, and so on. 
If the components under question are mass produced and their selection in a given system is assumed 
to be made at random, then the strength should also be considefed a random variable. The reliability 
of a component during a given period [0, T\ is taken to be the probability that its strength exceeds the 
maximum stress during the entire interval. The above formulation of the reliability problem was 
first proposed by Birnbaum [1], who also gave a method for estimating the related probability. In 
recent years, the model has found an increasing number of applications in many different areas, 
especially in the structural and aircraft industries. 

The stress to which a component is subject in its operating environment can often be predicted 
using the available technological knowledge about the mechanical, thermal and other relevant condi- 
tions of the system, and the way they interact with each other. In many cases mathematical models 
are available which allow one to compute deterministically the values of "stress" at different points of 
time corresponding to a given set of input values describing the initial conditions of the system. For 
example, as Church and Harris [2] have pointed out, in a missile flight these initial values may cor- 
respond to the propulsive force, angles of elevation, changes in atmospheric conditions, etc. Since 
these input values have statistical distributions associated with them, the calculated "stress" values 
will also be randomly distributed. Moreover, since it will be possible to generate a practically unlimited 
number of such stress values using simulation methods by varying the input values in accordance with 
their respective distributions and hence to reconstruct the stress distribution as accurately as desired, 
we might assume that the distribution of the maximum stress is known. The strength distribution can- 
not, however, usually be computed from such a priori considerations. It can only be estimated using 

159 



160 M. MAZUMDAR 

statistical methods from the results of tests carried out for this purpose. The number of such tests is 
likely to be small in most cases due to practical limitations. 

Under, the assumption that the distributions for strength and the maximum stress are both normal, 
we derive in this paper some point and interval estimates of reliability using "interference theory," 
when as above the stress distribution is known and the strength distribution is unknown, but a few 
observations on the latter are available. The point estimates are given in section 2, and the interval 
estimates are derived in section 3. 

Owen, Craswell, and Hanson [7] considered a similar problem and gave confidence intervals for 
the reliability when both stress and strength are (a) jointly bivariate normally distributed and observa- 
tions are in pairs, or (b) when they are independently normally distributed, but they have a common 
(but unknown) variance. Most recently, Church and Harris [2] considered the same problem and 
provided a procedure for obtaining confidence intervals, which they compare with that of Govindarajalu 
[3]. Their procedure, although empirically demonstrated to be superior to that of [3], is, however, 
inexact since it uses the asymptotic normal approximation of a given statistic and requires the substi- 
tution of the population mean and standard deviations by their observed sample values. Our procedure 
for obtaining the confidence interval is based on that given in [7] and, though exact in a certain fashion, 
is admittedly inefficient. Some comparisons of this procedure with that of Church and Harris [2] are 
given in section 3. Sampling experiments show that the average length of the interval associated with 
the Church and Harris procedure is smaller than that of our procedure although it appears that, in 
general, it fails to cover the true parameter value more often. 

2. POINT ESTIMATES 

Let X denote the maximum stress to which the component is subject over a given period of time. 
We assume that X is normally distributed with mean /jli and variance cr'f, and that both these values 
are known. Let Y denote its strength. We assume that Y is also normally distributed with unknown 
mean /x-z and variance erf. In what follows, we shall first assume that cr| is known and then that it is 
unknown. No doubt the latter is the more realistic case. In terms of our definitions we readily see the 
reliability to be 

(1) p = Pr{Y-X>0} 

= <b( ^-^ ) 

\ Vo-f + o-f/ ' 

where 

®(t)=—^= (' e- x *' 2 dx. 

Let Yi,Yz, . . . , Y m denote m mutually independent observations on the strength of the components. 
These values are obtained from bench tests. On the basis of these observations, we obtain below the 
minimum variance unbiased (MVU) and the maximum likelihood estimates of p. 

2.1 Case 1: cr| known 

Let Z be a normal random variable with mean and variance erf, and let 

(2) / = 



1 


if 


Y x + Z> /x, 





if 


Yi + Z^m 



RELIABILITY BY USING INTERFERENCE THEORY 161 

Then 

(3) E{I}=Pr{Y l +Z> f x 1 } 



=P- 

It therefore follows from Blackwell-Rao and Lehmann-Scheffe theorems (See Chap. 3. [5]) that 

_ _ »' 

the conditional expectation E{I\Y} is the MVU-estimate of p, where F=VF;/m is the sufficient 

i= 1 

statistic. Since moreover the conditional distribution of Y t + Z given F is normal with 

(4) £{F,+Z|F} = F, 

Var{F 1 +Z|F} = cr2 + o-|('l--) , 

the MVU-estimate of p in this case is given by 

(5) p, = E{h\Y) 

= Pr{Y i + Z>p. 1 \Y} 

= *( ^ ). 

VVo-f + o-fd-l/m)/ 

The maximum likelihood estimate of p, on the other hand, is given by 

It is seen that as m becomes large, the two estimates become equivalent. 
2.2 Case 2: erf unknown 

In this case the sufficient statistics is the pair (F, S.J) where 

»' _ 

Si=2 (F,-Y)*. 

Therefore, using arguments identical to those of Case 1, it follows that in this case, E{I/Y, S'i} is the 
MVU-estimate of p. From [6] we find the conditional density of Y\ given F and Sf to be 

' m \ 2 ) 1 f (Y-yV m ] 2 




m — ,- „/m — 2\ S-» 1 \ S-> J m — 



otherwise. 



162 M. MAZUMDAR 

Therefore the MVU-estimate of p in this case is 
(8) p n = E{I\Y,Si} 

= Pr{ Y x +Z >/xi | Y,S%} 



V m 



I i Y—y / m 
By use of the transformation 2 = - + - — - — \ - the integral on the right hand side of (8) reduces to 

z z 02 v m — I 



Jo r 



(9) />// = I ' . r( o m , 2) . zf"-'"/ 2 (l-z) (, "- 4)/2 ^(« + j82)ffe, 

'm — Z\ v (m — Z x 



where a =(s.,J— 1- Y—pn )/cr,, 



m 



and /3 : 



2S 2 \m — 



<j\ \ m 



It does not appear that this integral can be evaluated in a closed form; however, it can be evalu- 
ated very easily for particular examples by the use of numerical integration. 
Thf* maximum likelihood estimate forp in this case is 

(10) P// = <*>'- 



where 



Vaf+4 
S? 



ss = 



2 m-r 

pu is biased, but might be quite efficient when m is reasonably large. Also, it can be shown that for 
large values of m, the MVU-estimate pu and the maximum likelihood estimate/?// are equivalent. 

3. INTERVAL ESTIMATES 

As in the discussions in the preceding sections, we give the interval estimates of reliability for two 
separate cases (a) when cr\ is known, and (b) when o" 2 is unknown. We derive here only the lower 
confidence limit because of its importance in the present application. An extension of our procedure 
to two-sided intervals will be obvious. 

3.1 Case I: o~ 2 known 

Let K a be such that <$>(K a ) = a. Since Y is normal with mean p.-, and variance cr 2 /m, it follows that 



,11) Pr{ V "<r-*> «-*.}-!-„, 

f fii — fiz Ui — Y—K a <T 2 lVm} . 
or rr \ , =£ \ = 1 — a, 

I V o-? + o- 2 V rr? + 0-2 J 



'oi + 0- 2 



RELIABILITY BY USING INTERFERENCE THEORY 163 



Pr 



k <^1 h- 



Thus the 100 (1 — a) percent lower confidence limit for p in Case I is 

(12) p,,dl-a)=<t>l ./-—-; )• 

\ Vo-f + crj ' 

3.2 Case II: &\ unknown 

The procedure suggested by Church and Harris [2] for obtaining the interval estimate in this case 
is based on the statistic 

v- ? -^ 



^ai + sl 



They have shown that V is asymptotically normal with mean (ju, 2 — Pi )/ V^f^ and variance 

(/Lt 2 — (JLi) 2 <ri 



OV — 



m 2(m-l)(o-'f + o-2) 2 



1 ■ "2> 

Using normal approximation for the distribution of V and upon replacing the unknown values of fx 2 
and <t\ by their respective sample estimates Y ands*, they obtain (the symbol = denotes approximately 
equal to) 

(13) Pr { / *IZ^ > V+ K a CTv] = l-«, 

I Vo-2 + O"! J 

or Pr{p>4>(K+^„o-v-)} = 1-a, 

or Pr{p><D(r+^)}=l-a, 

where 



&l ai + ,\ 



■+ 



m 2(m-\)((T 2 + siy 



Since this procedure involves two approximations both of which may not be close in the non-asymptotic 
case, the above results may not be satisfactory in small samples. 

We next suggest a procedure based on results given in section 4 of [7]. This procedure is "exact," 
but is admittedly inefficient because it does not make use of all the available information; however, the 
use of sampling adopted in this method seems to be the only way we can eliminate the nuisance param- 
eter &\ without introducing any approximations, and we believe that the properties of the procedure 
might therefore be of interest. It seems possible that ad hoc procedures can be developed which use 
this sampling concept, but are more efficient than the procedure described here. Results obtained from 
simulation experiments given below seem to indicate that our procedure covers, in general, the true 
parameter value more often than the Church and Harris procedure, although the average interval 
obtained from the latter procedure is significantly shorter. 

According to our procedure, we first draw, using random number tables, a sample of m independent 
observations form a normal distribution with mean and variance cr\. Having taken these m observations 



164 M. MAZUMDAR 

we consider the transformed observations 



It is clear that 
(14) 

and 

Now let 

(15) 



and 

Then 
(16) 



Wi = Yi + Zi, i = l, 2, . . ., to. 



E{Wi}=p 2 



Vai{Wi}=(r\ + trl 



w=^ — 



TO 



£ {Wi-wy 



m^\ 



Sw 



Sw 



Vof + o-f Vo-, 2 + o- 2 2 J/ VqHd/ 



is distributed as a non-central < with degrees of freedom = to— 1 and the non-centrality parameter 
8— Vto(ju,i — /A2)/ Vo-f* + o"| . Let the left hand side of (16) be denoted by T,,, 1. It then follows from the 
formulas given in [4] that 

(17) Pr{ ^ (/"-**»> >S(m-l < r ffl - t , q) j = l-q, 

( Vo- 2 + crl J 

where from Eq. (33) of [4], 8(to— 1, T m -u ot) is defined as 

(18) 1 8(to — 1, r m _i, a) =T m -i — \(m— 1, r m _i, a) 



1 + 



71- , 



2(to-1) 



1/2 



and k(m— 1, r m _i, a) can be obtained from Table IV of Johnson and Welch [4]. From (17) and (18) 

it now follows that 

5(to — 1, T m -u a) 



(19) 



Pr 



{p>^( 



TO 



= l-a, 



so that a lower (1 — a) percent limit forp following this procedure is given by <J>(S(to— 1, T m - U a)/Vm). 



RELIABILITY BY USING INTERFERENCE THEORY 



165 



We now exhibit the results obtained from 200 artificial sampling experiments comparing the 
above two procedures using the same set of random numbers. In Table 1 the Church and Harris 
procedure is denoted as Procedure A and the method suggested in this paper as Procedure B. The 
particular sample sizes shown in the table were chosen because of the comparative ease with which 
the values of 8(m — 1, T m -i, a) can be evaluated from the existing tables. 

Table 1. Comparison of Confidence Interval Procedures 







0*i= 


0, o-, = l; 1 


-a = 0.95) 






p 


p-z 


0"2 


m 


Number of one-sided con- 
fidence intervals out of 200 
which covered true param- 
eter value 


Average length of 

200 confidence 

intervals 


Proc. A 


Proc. B 


Proc. A 


Proc. B 


0.90 


1.43287 


0.5 


10 


184 


191 


0.15 


0.28 








37 


183 


195 


0.12 


0.19 




1.81246 


1.0 


10 


171 


187 


0.19 


0.26 








37 


185 


193 


0.14 


0.19 




4.05278 


3.0 


10 


187 


183 


0.25 


0.27 








37 


187 


191 


0.17 


0.18 


0.95 


1.83905 


0.5 


10 


181 


188 


0.087 


0.19 








37 


188 


190 


0.067 


0.10 




2.32624 


1.0 


10 


184 


194 


0.12 


0.21 








37 


184 


190 


0.08 


0.11 




5.20163 


3.0 


10 


187 


189 


0.18 


0.20 








37 


185 


189 


0.10 


0.11 


0.99 


2.60055 


0.5 


10 


178 


187 


0.022 


0.095 








37 


186 


193 


0.015 


0.038 




3.28946 


1.0 


10 


182 


194 


0.045 


0.097 








37 


183 


189 


0.022 


0.039 




7.35546 


3.0 


10 


190 


193 


0.093 


0.11 








37 


192 


188 


0.034 


0.039 



REFERENCES 



[1] Birnbaum, Z. W., "On a Use of the Mann-Whitney Statistic," Proceedings of the Third Berkeley Symposium on Mathematical 

Statistics and Probability (Uni\ersity of California Press, Berkeley, Calif., 1956), Vol. 1, pp. 13-17. 
[2] Church, J. D. and Bernard Harris, "The Estimation of Reliability from Stress-Strength Relationships," Technometrics, 12, 

49-54 (1970). 
[3] Govindarajalu, Z., "Distribution-Free Confidence Bounds for P{X < Y}," Annals of the Institute of Statistical Mathematics, 

20, 229-238 (1968). 
[4] Johnson, N. L. and B. L. Welch, "Applications of the Non-Central t Distribution," Biometrika, Vol. .31, 362-389 (1940). 
[5] Lehmann, E., Notes on the Theory of Estimation, Mimeo Notes (University of California Press, Berkeley, Calif.) (Reprinted 

1962). 
[6] Lieberman, C. J. and G. J. Resnikoff, "Sampling Plans for Inspection by Variables," Journal of the American Statistical 

Association, 50, 457-516 (1955). 
[7] Owen, D. B., K. J. Craswell, and D. L. Hanson, "Non-Parametric Upper Confidence Bounds for P{Y < X) and Confidence 

Limits for P{Y <X} When X and Y Are Normal," Journal of the American Statistical Association, 59, 906-924 (1964). 



THE OPTIMAL BURN-IN TESTING OF REPAIRABLE EQUIPMENT* 



John M. Cozzolino 

Department of Statistics and Operations Research 
University of Pennsylvania 
Philadelphia, Pennsylvania 



ABSTRACT 

The subject of this paper is the utilization of the "infant mortality" or decreasing failure 
rate effect to improve the reliability of repairable devices. Decreasing failure rate implies 
the possibility that devices which exhibit it can be improved by "burn-in testing" of each unit. 
Such a test serves to accumulate operating time while shielded from the full costs and con- 
sequences of failure. A general formulation of the burn-in test decision for repairable devices 
is presented and some special cases are solved. A class of models, indexed by the degree of 
partial replacement present in the repair process, is considered and numerical results for 
the optimal policy are given for several members of that class. A comparison of those results 
reveals the profitability of testing increases with the complexity of the repairable device. 



1. INTRODUCTION 

The subject of this paper is the utilization of the "infant mortality" effect to improve the reliability 
of repairable devices. The term infant mortality denotes a decrease with age of the conditional proba- 
bility of failure of a device and, consequently, an improvement in reliability with time as the device 
operates. This phenomenon is widely observed. For example, decreasing failure rate has been found 
to occur in electronic systems and components [12], biological systems [5], social systems such as 
business enterprises [20], and also, in mechanical systems. In the case of repairable equipments de- 
creasing failure rate is often observed after the occurrence of failure and repair. 

Decreasing failure rate implies the possibility that devices which exhibit it can be improved by 
"burn-in testing" of each unit. The term burn-in testing is taken to mean test stand operation under 
environmental conditions as close as possible to true operating conditions. The basic purpose of this 
test is to accumulate operating time while shielded from the full consequences and costs of failure. 
The measurement of reliability, while not a primary goal, may be a partial by-product. 

The main theme of this paper is a general formulation of the burn-in test decision problem for 
repairable equipment. A continuous stage dynamic program is used to treat this sequential decision 
problem and some solutions are given for special cases. Space does not permit a full discussion of the 
class of models for which this decision structure applies; however, one broad family of models, indexed 
by the degree of partial replacement present in the repair process, is presented. Numerical results 
for the optimal test policy are given for several members of that family, and a comparison of those 
results reveals an increase in the profitability of testing as a function of the complexity of the repairable 
device. 






*The work herein reported was supported by the Office of Naval Research under Contract Nonr 1841 (87); NR 042-230 while 
the author was at the Operations Research Center at the Massachusetts Institute of Technology. The computation reported here 
was done in part at the Computation Center at Massachusetts Institute of Technology. 

167 



168 J. M. COZZOLINO 

2. LITERATURE REVIEW 

The distribution theoretic conditions under which the expected future life of a device increases 
with its age are shown by Watson and Wells [27]. Based upon the work of Knight [18], one stated 
result is that the expected future life of an aged device always exceeds that of a new one if the device 
has a monotonically decreasing failure rate. 

Several authors have investigated the burn-in test problem for the case of unrepairable equipments. 
In that problem a failed device is considered worthless. The quantities investigated are the burn-in 
time required to achieve a specified expected future life of the surviving units. Watson and Wells [27] 
discuss these quantities for several failure distributions. A nonparametric method is considered in 
a paper by Lawrence [19]. Assuming only a monotonically decreasing failure rate, he derives bounds 
upon the burn-in time required to achieve a specified expected future life. 

Recently, interest has developed in the burn-in testing of repairable devices. This potentially 
important problem requires some explicit assumptions about the nature of the relevant failure and 
repair processes. This necessity follows from the need to assess the relative worth of the repaired 
device. Kooharian and Young [28] consider a useful class of reliability models for repairable complex 
devices. The outstanding feature which they prove for this class is a special simplified form of the 
burn-in test decision policy. A similar model for the "debugging period" is given by Bazovsky [4]. 

Concern for the long run reliability improvement of complex systems is evidenced by a paper by 
Barlow and Scheuer [3] entitled "Reliability Growth During a Development and Testing Program." 
They attributed the improvement mainly to engineering changes. 

In a report upon which this paper is based [7], the author presents some models of decreasing 
failure rate processes incorporating explicit consideration of the repair process. This union of failure 
and repair models appears to be a fruitful basis for the discussion of burn-in testing in view of the 
interaction between the two processes. A complete picture of the life of the device is given in terms 
of a stochastic process, usually with non-independent interfailure times. 

3. ELEMENTS OF THE PROBLEM 
The Unit Test 

A unit test starts with the testing of a new unit and ends with the output of one tested unit. If a 
failure occurs during test, the test is interrupted, the unit is repaired, and the test then resumes. The 
time elapsed during repair is not counted. 

The Failure Process for New Units 

Given a population of indistinguishable new devices, a specified and constant environment in 
which they are intended to operate, and a set of conditions which define failure, a distribution function 
Fo(tl6) with parameters 6 and density fo(t/6) is utilized t6 represent the uncertainty which exists re- 
garding t, the duration of operating time prior to the occurrence of failure for any one of the devices. 
The failure rate of that density, also known as the hazard rate, conditional probability of instantaneous 
failure, or force of mortality, is a very important function for determining the reliability characteristics 
of the device. For a distribution satisfying Fo(0~/6) —0, the failure rate, a function of time, is defined as 

It will be assumed throughout thatfo(t/6) is a continuous function of t and that Fo(O~/0) = 0. 



OPTIMAL BURN-IN TESTING « 169 

Density of Future Time 

The original density, /o, represents the process only before testing begins. Suppose that at some 
time, r, after the test starts no failure has been observed in (0, t). Let t = future time to failure. Then 
the total time to failure is (£+ r). Of this total time to failure the only random part remaining is t. The 
relevant density function is that of future time to failure. This renormalized form of the density will be 
denoted here by 

f mr)- M,+Tle) - 



[l-F„(r/«)] 

While this density function usually represents a larger family of densities than does/oU/0), its failure 
rate is the same except displaced to the left by t; no(t + Tl8). 

The Repair Process 

We consider that a failure during test will be followed by a repair of the failed unit. The rule which 
assigns a new density of time to next failure for the repaired unit will be called the "repair rule." 

In many analysis it is assumed that the repaired unit has the same density of time to next failure 
as the new units. This practice ignores the information available from the failure and from the past 
history of the unit. Since it is possible to distinguish between new and repaired units we realize that 
it may be valuable to do so. This is particularly relevant for burn-in testing models since the purpose 
of testing is to create a difference between new and aged units. 

The only restriction which we shall place upon the new density is that it be from the same family 
of density functions as is f{t/8, t). Hence we allow t to have the meaning of a parameter of the density 
rather than strictly the time since last failure. 

The repair rule then is a transformation which determines which density of that family applies. 
The most general transformation to be considered here is one that allows the new parameter values 
to be a function of the parameter values which applied just prior to the failure. Hence the new parameter 
values are 

d R =o R (o,T) 

and 

tr = t h {6, t), 

where (8, r) are the parameter values just prior to failure. 

These functions will be required to have the property that, when applied to any point in the space 
of allowed parameter values, they will produce a point which is also inside that space. Also, their con- 
tinuity in t will be required. 

The repair rule results from an explicit model of the repair process. Examples of this are given 
in [71 

The Reliability Value Function 

The reliability value function specifies the value that the tested unit would have if testing were 
terminated. It gives the value as a function of any possible test position (6, t). It will be denoted by: 
U(8, t) = the values of a unit when testing is terminated and a time when the state of process is {8, t). 

This class of functions, U(8, r), includes all expectations of the future behavior of the unit, since 
it includes all the parameters of the density function of future life. 



170 J- M. COZZOLINO 

The Cost of Testing 

The last remaining element of the problem is the cost of testing. In this we seek to represent 
the actual cost of using the test facilities. The cost of testing will be considered to consist of the fol- 
lowing two components: 

Ci = The cost of repair of a unit which fails during test. When several independent failure modes are 

considered this definition can be expanded to include a different cost for each failure type. 

C2 — The cost per unit test time. 

4. OPTIMIZATION BY DYNAMIC PROGRAMMING 

The profit from a unit test is the terminal value (given by the reliability value function) minus 
the test costs incurred during the testing process. At an intermediate time during the test the optimal 
test policy does not depend upon the costs already incurred since they are fixed or "sunk" costs. 
The optimal policy is that which maximizes the expected profit of the remaining portion of the unit test. 

The decision-making process is sequential; the decision to stop or continue is remade after each 
observation of the partial test results. Since partial test results are revealed continuously in time each 
time point is, in fact, a decision point. The basic repeated unit of the decision tree is shown in Figure 1. 



/ / 



, DEC IS ION NODE 
STOP TEST FAILURE 



NODE 



-•-O CHANCE NODE 
NO FAILURE 



DECISION NODE 

FIGURE 1. The basic unit of the Decision Tree 
We define the optimal value function as 

V(6, t) =the expected value of the remaining portion of the unit test when the current state is (6, t) 
and when all future test decisions are optimal. 

Let us consider a small positive time interval, h, starting at a time when the state of the process is 
(6, t). The optimal value function, V(6, t), can be expressed in terms of the optimal decision at the 
start of the interval, the events which can occur during it, and the optimal value function at the possible 
resulting states. Let us consider the possibility of one event occurring in (t, t + A) at time r + e, where 
e < h. 

Since there are two possible decisions the optimal' value function is equal to the maximum of 
two terms. The equation is 



V(0, T) = Max{tf(0, t),-C 2 A + 



F (T + hlO)-Fo(Tldy 



1-Fo(t/0) 



[-d + VlOniO, r+e), t h (6, r+e))] 
l-F (T+h/0) 



+ 



l-Fo(rld) 



V(d,T+h) 



The term V(d, t) can be subtracted from both terms in the Max bracket. The resulting expression can 
be written as 



0=Max (J(0, t)-V(0, t), h \-C 2 



,t),a{-< 



OPTIMAL BURN-IN TESTING 171 

V(d, T+h)~V(d,T) 



+ 



F (T+hld)~Fo(Tld) 



h(l-F o (Tl0)) 



[-C l + V(6 R {0,T + e),T R {d,T + e))-V(0,T+h)]\\. 



Let us now observe that we have the form 

O=Max{£/(0, t)-V(6, t), h£(h)}, 

where £(h) represents the interior of the second term. It should be noticed that the value of A has 
no effect through the term £(h) since when hi;{h) <0 the function V(6, t) is determined by the 
equation 0= U(0, t) — V(6, t). In fact, for any positive h, 

h£(h)<0 <^f(A)<0,Af(A)=0<^f(A)=0. 

Therefore, we can remove h and write the equation as 

O = Max{£/(0, t)-V{6, t), £(h)}. 

The interval, h, can be considered as an inspection interval with the understanding that if a failure 
occurs anytime during the interval, its occurrence will not be detected until time, t+ h, at which time 
the operation of the device will be resumed. If detection is not an important aspect of the problem, 
however, the limiting form of the equation is more useful. Provided only that the functions /o(t/0), 
6h(6, t), and tr{6, t) are continuous functions of t everywhere in their range of definition, we can 
now let h approach 0. In the limit, we obtain the equation 

0= Maxjf/(0, t) - V(6, r) , dV{ ^ T) - «o(t/0) [V(6, t) - V(d«, r«) + C,] -cX 

It is helpful to visualize the testing process as a point moving with time in the (0, t) space. As 
testing proceeds, r (effective age of the unit) increases and remains constant at non-failure points 
in time. When a failure occurs, both and r normally jump to a new value given by (0«, Tr). The 
solution of the dynamic program divides the (0, t) space into two types of regions: 

1. Continuation region where 
(CI) U(6, t)-V(6, t) «0, and 

(C2) ^ ( ^ t) -/xo(t/0)[F(0, T)-r(0 ft ,T H ) + C,3-C 2 = O. 

2. Termination region where 
(21) U(d, r)-V(d, t)=0, and 

(T2) d ^ T) - M o(T/0) [^(0, r) - F(0«, r«) + C,] - C 2 ^ 0. 

These regions are separated by a boundary, t*(0), which represents the decision to terminate 
testing. A simple illustration of a process trajectory in state space is given in Figure 2. 



172 



J. M. COZZOLINO 



STATE SPACE (ONE 9 PARAMETER) 

•(e) 



REPAIR JUMP 




TERMINATION 
BOUNDARY 



o Termination Point 
X Failure Occupied 



FIGURE 2. A process trajectory in state space 

The following theorem helps to characterize the boundary in terms of the properties of U(6, r). 
THEOREM: For each 6, the optimal policy, t*(0), is restricted to the following possibilities: 

a) t* = or oo. 

b) t* = t such that the derivative, 2 , is discontinuous at t. 

St 

, * . , dU(6, t) . 

c) t =t such that 2 is continuous at t, 

dr 

d U(0, t) 8V(6,t) , a 2 t/(6>, t) d 2 V{d, t) 

, and — s= — at r. 



8t 



dr 



dr 2 



3t 2 



PROOF: Assume that, for given 6, V{6, t) and t*(^) are the optimal value function and decision 
boundary, and that 

a) 0<t*(6>) <oo 

b) U(6, t) has continuous first derivative, - 1 at t*. 

OT 

At a point {6, r* — e) (for e > 0) in the neighborhood to the left of t*, both V(d, t) and U{6, r) and 
their derivatives exist by assumption. Condition CI gives 

U(6,T*-e)-V(d,T*-e) ^0. 

Expansion about the point t = t* gives 



U(6, T*-e) - V(6, T*-e)=U(d, r*) - V(d, r*) -e 



dU(d,r*-ye) dV(6,T*-ye) 
dr St 



for some < y < 1. Conditions CI and T\ then imply the following condition in terms of the left-hand 
derivatives 

dV-{e,T*) ^ dU(e,T*) 

dr dr 

The left-hand derivative also satisfies Eq. C2 

ay- (e, t*) 



dr 



fjLo(T*/e)[V(e, T *) - v(o R , t«) +c] -c 2 =o. 



In a neighborhood to the right of r* the conditions 71 and T2 give 

dU+(d, T*) 



dl 



■ti (T*id) [U(e, T ) - v(o K , th) +C] -c 2 ^ 0. 



OPTIMAL BURN-IN TESTING 
Subtracting the latter two conditions and using 71, we have 

dV-(e,T) dU+(e, t) 



173 



dr 



dT 



Assumption C ensures that both one-sided derivatives of U are equal; hence we have the existence 
of Ps derivative at t* and the equality 



dV(e, r) dU(e, t) 



dT 



dr 



atr = r. 



Furthermore, if the functions U and V have second derivatives at t* then a Taylor series expansion 
using both derivatives and the above equation yields 



d 2 V(d, r) d 2 U(d, t) 



dT 2 



dT 2 



at r = r*. 



Hence part Q of the theorem is proved. The possibilities excluded in assumptions a and b are parts 
a and b of the theorem. 

The differential equation is a necessary condition upon the optimal value function. Its form and 
complexity will depend upon the nature of the 6 space and repair functions. In the simplest cases the 
complete family of solutions of the differential equation can be obtained analytically. The parameters 



V(8,r) for fixed 8, 


| 




1 

1 ___ 

L_- ---" V(8,t) = U(8,t-) 


V(8,r> _______ ; - ; ^^= , •' 


1 


_____ U(8,t) 


1 
| 


CONTINUATION REGION 


I TERMINATION REGION 



r*(8) 
SOLUTION AT TANGENCY 



V(8,r)FOR FIXED 8 




CONTINUATION 
REGION 


| TERMINATION REGION 

1 

V = U ________ 


V(8,r> ^_____- :5 ==>-^' 


1 


^^^Ui6,r) 


1 

1 
1 



SOLUTION AT DISCONTINUITY 


OF DERIVATIVE 


V(8,t)F0R FIXED 8 






CONTINUATION REGION 




TERMINATION REGION 


V(8,r) ^ — - 




V=U 




U(8,r)-0 







SOLUTION AT DISCONTINUITY OF U(8,r> 

Figure 3. Possible types of solution 



174 J. M. COZZOLINO 

of its solution can then be evaluated using the above boundary conditions. This will be illustrated later 
for two special cases. Figure 3 illustrates the tangency condition and the two types of discontinuity 
conditions. 

The method of optimization discussed here is a continuous stage dynamic program. Other mathe- 
matically similar problems in the literature, such as that of age replacement policy for stochastically 
failing devices which wear out, have been formulated as a dynamic program using a "renewal equation" 
form. For example, [11], [15], [16], and [17] use that form. Either form could be used here although this 
differential form is thought to have some computational advantage for this problem. 

The Distribution of Remaining Test Time 

When analyzing the feasibility of burn-in testing it is useful to know certain statistics of the unit 
test. Perhaps the most important of these are the moments of the test duration which will be con- 
sidered here. Let t be the remaining test time when the test process is in state (0, t), D(t/6, t) its 
density, and D(S/6, t) be the Laplace transform of the density. An equation for D{Sld, t), similar 
to the value function equation, can be derived using similar argument. By considering a small time 
interval, h, and the events possible within it, the following equation is derived: 

D(Sld, t) = e-™ D(Sld H , t r ) ixti(rld)h + e- sh D(SI6, t) [1 -/lh>(t/0)A] for t<t*(9). 

By allowing h to approach zero, the limiting form obtained is 

dD(Sie,T) _^ {T/e) [D(Sldr) _ DiS i eRTH)]=SDiS/dT) for T<T * ( 0). 

The associated boundary condition is 

D(S/6,t*) = 1. 
As a consequence, the Arth moment of remaining test time, defined as 



i 



E k (d,T) = \ t k D(t/d, r)dt for it = 0, 1, 2, . . ., 



obeys the following equation: 
dE k (0, r) 



-/ao(t/0) [E k (d, t) -E k (0 H ,T R )] =-kE k .t(8, t) 



dT 

for A=l, 2, . . ., and t 1 t*(0). The boundary condition is 

E k {6, t*)=0 for A=l, 2, 3 

The solution for these moments requires the knowledge of t*(0); then, since the (k+\) st equation 
involves the Ath moment, the equations for the moments can be solved successively starting with 
A=l. 

5. EXTENSIONS OF THE FORMULATION 
The Optimal Policy Under Discounting 

In the analysis of processes which take place over an appreciably long time, thus tying up capital 
equipment and possibly delaying cash flows, it may be useful to account explicitly for "value of time" 
by optimizing the present discounted value of the process. Under this method, the values of all future 



OPTIMAL BURN-IN TESTING ' 175 

value transactions are discounted from face value according to the length of time into the future at 
which they take place. Hence, if the amount, £/, is to be received at time, T, into the future, its present 
discounted value is Ue~®'\ where /3 is a positive discount factor. 

The optimal policy and value function under discounting will be different from that already dis- 
cussed; however, only a slight modification of the optimization equation is required to treat this case. 
The equation now appears as 

0= Max {U(d, t) - V(6, T), dV( ° ,T) -pV(6, r) - u (tI6) [V(0, t) - V(0 H , r^-C,] + C 2 }- 

OT 

The one additional term is (3V(6, t). All that has been said about the non-discounted optimal solution 
applies to this case also. 

The effect of present value discounting will usually be to shorten the test. This can be expected 
since the terminal value is always received at the end of the test, whereas the test costs are spread 
out over the duration of the test. Therefore the expected present discounted value of a given policy is 
necessarily less than the expected undiscounted value of the same policy. 

Burn-In Testing Under Uncertainty of the Failure Process 

We have so far considered a known failure process. However, in many applications there exists 
a high degree of uncertainty about the underlying failure process. One method of relaxation of this 
assumption of complete knowledge of the density is to consider a known form of density function but 
with some unknown parameters. The state of knowledge about the values of such parameters is ex- 
pressed, at any given time, by a joint density function over all allowable sets of their values. This prior 
density is modified by the observation of results of the testing process according to Bayes' Rule. The 
concept of a conjugate prior density, as described by Raiffa and Schlaiffer [25], is used here to enable 
a mathematically tractable formulation of this problem. 

We suppose that the parameterization of the density includes parameter(s) cf> which are unknown, 
but do not change with failure. They are now included explicitly in the density, f(t/6, t, </)). The "prior" 
density, that which describes the state of knowledge at a given time, is denoted g((/)/i//) . It is normalized 
over the range of possible <f) values for any values of its parameters, t|/. 

The density of fife, unconditional upon the unknown 0, is 



Jail <* 



h(tie, r, <//) = f(t/e, r, 4>)g{<t>l<\>)d<i>. 

Jail <b 

The prior density is chosen to be conjugate to f(t/<f), t, </>) so that, after observation of a failure at 
time, t , following start-up of a unit of effective age, t, the posterior density, according to Bayes' Rule. 

/(tie, r, c/QgW) 

h(tld, t, i//) 

has the same functional form in as does g{<j), <//), but with parameters changed to «//f = t///. (6, t, (//, t). 
Hence the value of the parameters of the prior density is sufficient to describe the state of knowledge 
at any time. 

The failure rate function which appears in the equation is 

V(tld, T, (//)= 



/; 



h(x/6, t, ty)dx 



176 J- M. COZZOLINO 

The state of the process under uncertainty is given by (6, t, i//, x) , where ^ andT are the values current 
at the time of last failure or start-up and x represents time since that event. Let this description be 
abbreviated by the definition 

d' = {6, t, ¥}. 
The state after failure in state (8' , x) is 

d' R {6\ x) = {6, t + x), Tr{6, t + x), V F {d, T, ¥, x)} 

x R (d', x) = 0. 

It is now apparent that the adaptive burn-in test problem, for one unit test, has the same form of opti- 
mization equation as previously described: 

0=Maxj £/(0\ x)-V(0', x), dV{ ®' ,x) -v(xld')[V(0', x)-V(d' R , 0) -Ci]+&] . 

It will, however, be more difficult to solve because of the larger state space. 

6. SOME SPECIAL CASES 

There are some important special cases of burn-in test processes for which there are no 6 type 
variables present; i.e., the failure density of a repaired unit is of the same form as that of some aged 
unit which never failed. In such cases, the effective age of the unit (t) completely describes the satte 
of the process and the differential equation holds on a line. We now consider three such cases which 
encompass a range of possible repair rule cases. 

Complete Replacement or Unrepairable Equipment (t« = 0) 

In the case of unrepairable equipment, a failure necessitates a complete replacement of tlv* .m J 
device by a new, untested one. Therefore the state of the process after failure is independent of the 
age attained by the old unit just prior to failure. This case is represented by the repair rule Tr(t) = 0. 

The resulting linear first order differential equation, 

^^-- (jlo(tI6)[V(t) -V(0)]^ fil{Tl6)C 1 + C 2 
is easily solved to give a family of all solutions parametric on one constant. For convenience, let 



w(t) =j- „ , -■ + 



C 2 f [l-Fo(xld)]dx 

Jo 



[l-Fo(Tld)] [l-F (Tld)] 

Since the constant of solution is determined by the boundary conditions we write it as t*. The solution 
is then 

( V{t*) - W(t*) + W{r) for ^ r ^ t* 
V(t)=\ 

[U(t) fort*<r. 

The value of t* is then determined as indicated by the previous theorem. 

The independence of interfailure times allows easy solution for the Laplace transform of the density 
of remaining test time. It is 

D{S,t) = Q(s',t*) ''*'*'* f ° r T ^ T * 



OPTIMAL BURN-IN TESTING J 77 

where 

1-J e- Sx Mx/e)dx 

Q^,r)=— [ [_ Fo(Tle)] 

Although the solution for this case could easily be obtained by more direct arguments, it is pre- 
sented here both to illustrate the method and to enable a comparison between unrepairable and re- 
pairable devices. We can already see that burn-in testing is somewhat inefficient for unrepairable 
devices since a failure results in the loss of all test effort expended up to the time of failure. 

Superficial Failure (tr — t) 

In the previous case, the test time "invested" in the current unit was completely lost with a failure. 
It represents one extreme case the opposite of which is called here superficial failure and results 
in no loss of invested time. This case, represented by the repair rule t«(t) =t, is important because of 
the simplification of the decision structure. 

The differential equation for this case is 

The optimal value function in terms of t*, is 

\U{t*)-C, \ fXo(x/d)dx-C 2 (T-T) forr^r* 
V(t)=\ h 

[U(t) for t s*t*. 

The Laplace Transform for the density of remaining test time is easily found to be 

Z)(S/T)=e- s(T *-"forTS=T*. 

This is the transform of a unit impulse function located at t = t; consequently the length of remaining 
test time is not random. The sequential decision which remains unchanged in the face of partial test 
experience. 

At first glance it seems counterintuitive that a failure leaves the device no better or worse than if 
no failure occurred; however, there are reasonable assumptions which lead to this repair rule. 
Kooharian and Young [28] present a class of such models. Similarly, Cozzolino [7] discusses a special 
case of that class, based upon the concept of production defects initially present in the device. Bazovsky 
[4] gives a similar debugging model. Because of its special deterministic nature, the decision policy 
is easy to implement. Also, as we shall see, this model results in a relatively profitable burn-in test. 

Partial Renewal (tr = (xt). 

The repair rule tr(t) =oct for 0<a<l represents an intermediate degree of replacement. 
Consider, for example, a device composed of N indistinguishable and essential components. Failure 
results from failure of any one of the components. Repair is then accomplished by replacement of that 
failed component by a new, untested component. The test time previously accumulated on all other 
components is retained. This situation may be represented by a partial loss of effective age, where 

a in the repair rule is a fraction of the time retained; I J in this case. A more detailed discus- 

sion of this model is given in [7]. 



178 



J. M. COZZOLINO 



The differential equation for this class of repair rules is 
dV(i 



-£- L -VL (Tl8)[V(T)-V(aT)]=Vo(Tld)C 1 +C 2 . 



The combination of differential and difference operators present makes analytic solution difficult for 
most failure rate functions of interest. However, approximate solutions based upon numerical integra- 
tion are easily obtained. Since all solutions are within an additive constant of each other, the boundary 
conditions are easily satisfied. 

A Specific Example 

Some numerical results, based upon the above family of repair rules, were obtained in order to 
illustrate how the optimal test policy depends upon the repair rule parameter, a. The specific model 
used is based upon the density function, 



| b t a \ 6+1 

Mtla, b)=- -3- *2=0, 6>0, a>0. 



The failure rate associated with this density is monotonically decreasing, 

b 



/xo(tla, b) = 



a + t 



The value of a tested unit is taken to be a function, R(t), of its first life after test, t, as shown in Figure 4. 



I 15 Ro' 

Ro' 

85 Ro 



RO) ' VALUE RECEIVED FOR TESTED UNIT 




NORMAL 


INCENTIVE RANGE 


PENALTY 


RANGE 




RANGE 







Figure 4. Value vs performance 



The reliability value function is consequently the expected value of this function as of test termination 
and is a continuous function of age, 



U(t) 



-\: 



R(t)f(tl6,r)dt 



= R« 



0.85 + 0.15 



a + T 



a + T + 7 7 , 



+ 0.15 



a + r 



a + T + T 2 



The parameter values chosen for the example are a = 72.0, 6=1.1, Ci = 1.0, C 2 = 0.01, R o = 90.0, 
ri = 500, and T 2 = 1000. The mean life is found to be 720 units and the value of an untested unit is 
U(0) = 78.57. The optimal test policy and value function were found for several values of repair rule 
parameter, a. The optimal policy, r*, is tabulated for these a values in Table 1, along with other quan- 
tities of interest. Figure 5 gives the optimal value function for one a value; a = 7/8. Figure 6 shows how 
the expected test profit, V(0) —U{0) increases with a. 



OPTIMAL BURN-IN TESTING 



179 



89 


- 




87 


- 


^ 


Sale Price 


85 








83 






>Test Cost 


81 








79 


/ 


■ 


>Tesl Profit 


77 


I 1 1 1 1 1 1 1 


1 1 1 1 1 1 



40 80 120 160 200 240 280 320 360 400 440 480 

FIGURE 5. The optimal value function V{t) for a = 7/8 



V(0)-U(0) 



18 






16 






14 






12 






1.0 






08 






06 






04 






02 






n 


- t — i — i— -~~\ i i i i i 




c 


) 01 02 03 04 05 06 0.7 0.8 09 1 


D 




6 8 10 = 


l-a 


i iii mi i 



Figure 6. Test profit as a function of a; 
Case 3 



The parameter, a, indexes a family of partial renewal models; the fraction of the device which is 
renewed with repair of a failure is (1 — a)= 1/iV. As a ranges from to 1 the model ranges from com- 
plete replacement (a=0) to the superficial failure case a— 1, and is therefore capable of representing 
a wide range of repair behavior. The computations show that as the degree of partial renewal decreases 
{a increases), the coefficient of variability of the total test duration decreases, as would be expected 
from previous results. As the test duration becomes more deterministic the optimal test policy, t*, 
lengthens. The expected profit due to testing increases and, in addition, the ratio of expected test 
profit to expected test cost increases, showing an inherently greater profitability as the renewed portion 
of the device decreases. This result conforms to our intuitive expectation based upon the previous 
results; as a increases a larger portion of the invested test time is retained at failure. Consequently, 
failure becomes less of a risk and the test process becomes closer to a deterministic one. 



180 



J. M. COZZOLINO 
Table 1. Optimal Test Policy for a Class of Repair Models 



a 


N=a-ar> 


Optimal 
policy 

T* 


Expected value 


Expected 

profit from 

test 

[V(0)-U(0)] 


Ratio (%) of 

test profit 

to test cost 

r V(0)-U(0) i 

IV(T*)-V{0)\ 


Expected 

duration of 

unit test 


Coefficient 

of variability 

of duration of 

unit test 


Before test 

V{0) 


After test 

V(T*) 





1 


14.5 


78.595 


78.979 


0.0249 


6.5 


15.7 


0.28 


0.5 


2 


80.0 


78.70 


80.68 


0.13 


6.6 


97.5 


0.26 


0.6 


2.5 


140.0 


78.83 


82.04 


0.26 


8.1 


174.2 


0.25 


0.7 


3.3 


220.0 


79.09 


83.63 


0.52 


11.5 


268.4 


0.20 


0.75 


4 


261.0 


79.27 


84.36 


0.70 


13.8 


310.9 


0.17 


0.8 


5 


304.2 


79.47 


85.07 


0.90 


16.1 


351.7 


0.14 


0.875 


8 


370.0 


79.82 


86.05 


1.25 


20.1 


406.4 


0.09 


0.9 


10 


390.0 


79.94 


86.32 


1.37 


21.5 


420.6 


0.07 


1.0 


oc 


470.3 


80.43 


87.36 


1.86 


26.9 


470.3 


0.0 



7. CONCLUSIONS 

The analysis of the burn-in testing of repairable devices requires explicit information about the 
repair process in order to evaluate repaired devices. A class of repair models, having a form tractable 
with respect to the sequential decision problem for optimization of the burn-in test, is described by a 
repair rule of the form, 0«(0, t), tr(0, t), specifying the parameters of a density function within the 
family of the renormalized form of the new unit density. This family of densities has the same func- 
tional form of failure rate. 

The burn-in testing optimization, when treated as a continuous stage dynamic program, results 
in a differential equation similar to equations for the state probability of a stochastic process. In two 
important special cases, the differential equation is of the first order in total derivatives and conse- 
quently is easily solved and fit to the boundary conditions. In more complex cases it has a differential- 
difference nature. 

The wide differences in the expected test profit under different repair models of the partial renewal 
family leads to the conclusion that the explicit modeling of the failure-repair process is highly important 
for the burn-in testing of repairable equipment. While this result is based upon a specific class of 
models, it appears likely that the greater profitability of burn-in testing of repairable over unrepairable 
devices may prevail over a much wider class of models. 

REFERENCES 



[1] Barlow, R. E., A. W. Marshall and F. Proschan, "Properties of Probability Distributions with Monotone Hazard Rate," Annals 

of Math. Stat., 34, 375-389 (1963). 
[2] Barlow, R. E. and F. Proschan, Mathematical Theory of Reliability (Wiley, New York, 1965). 



OPTIMAL BURN-IN TESTING - 181 

[3| Barlow, R. E. and E. M. Scheuer, "Reliability Growth During a Development Testing Program," Rand Memorandum 

RM-4317-NASA (Nov. 1964). 
[4] Bazovsky. L, Reliability Theory and Practice (Prentice Hall Inc., New Jersey, 1961). 

[5] Chemical Rubber Publishing Company, Standard Mathematical Tables, 12th edition. Commissioner's 1941 Standard Ordi- 
nary Mortality Table. 
[6] Cox, D. R.. Renewal Theory (Methuen and Company, Ltd.. London, 1962). 

[7] Cozzolino. J. M., "The Optimal Burn-In Testing of Repairable Equipment." Massachusetts Institute of Technology, Opera- 
tions Research Center Technical Report No. 23 (Oct. 1966). 
[8] Davis, D. J., "The Analysis of Some Failure Rate Data," J. Am. Stat. Assoc, 47, 113-150 (1952). 
[9] Drenick, R. E., "The Failure Law of Complex Equipment," J. Soc. Indust. Appl. Math., Vol. 8, No. 4 (Dec. 1960). 
[10] Dreyfus, S.. "Dynamic Programming and the Calculus of Variations," Rand Dept. R-441-PR. 
[11] Fox, B. L., "An Adaptive Age Replacement Policy." Operations Research Center, University of California, Berkeley 65-17 

(PR) (1965). 
[12] General Electric Co., Transistor Manual, Chap. 18, "Semiconductor Reliability," 6th Ed. (1962). 
[13] Haddon, M. O, H. Cheesebrough, and R. E. Kirby, "Guaranteed Reliability at a Minimum Cost," IEEE Trans, on Reliability 

and Quality Control (Sept. 1958). 
[14] Hausman, W. H., and M. Kamins, "Early Failures in Automobile Parts: A Background Study in Reliability," Rand Mem. 

RM-4002-PR (1964). 
[15] Jorgenson, D. W. and J. J. McCall, "Optimum Scheduling of Replacement and Inspection," Operations Research II, 732-746 

(Oct. 1963). 
[16] Jorgenson, D. W., and R. Radner, "Optimum Replacement and Inspection of Stochastically Failing Equipment," in Studies 
in Applied Probability and Management Science (Stanford University Press, Sanford, 1962), edited by K. J. Arrow, S. Karlin, 
and H. Scarf. 
[17] Kamins, M. and J. J. McCall, "Rules for Planned Replacement of Aircraft and Missile Parts," Rand Mem. RM-2810-PR 

(Oct. 1961). 
[18] Knight, W. R., "Exponential and Subsexponential Distributions in Statistical Life Testing," Ph. D. Thesis, University of 

Toronto (1958). 
[19] Lawrence, M. J., "An Investigation of the Burn-In Problem," Technometrics, 8, 61-71 (Feb. 1966). 

[20] Lomax, K. S., "Business Failures: Another Example of the Analysis of Failure Data," J. Am. Stat. Assoc, 49, 847-852 (1954). 
[21] Markowitz, D. H. and R. G. Pouliot, "Infant Mortality Screening; A Discussion of Burn-In and Other Screening Techniques," 

U.S. Naval Weapons Quality Assurance Office, Washington, D.C. (June 1967). 
[22] Morse, P. M., Queues, Inventory, and Maintenance (Wiley, New York, 1958). 
[23] Norris. R. H., "Run-In or 'Burn-In' of Electronic Parts," Proceedings of the Ninth National Symposium on Reliability and 

Quality Control (1963), pp. 335-357. 
[24] Proschan, F., "Theoretical Explanation of Observed Decreasing Failure Rate," Technometrics, Vol. 5, No. 3 (Aug. 1963). 
[25] Raiffa, H. and R. Schlaifer, Applied Statistical Decision Theory (Harvard University Press, Cambridge, 1961). 
[26] Riordan, J., Stochastic Service Systems (Wiley, New York, 1962). 
[27] Watson, G. S.. and W. T. Wells, "On the Possibility of Improving the Mean Useful Life of Items by Eliminating those with 

Short Lives," Technometrics, 3, 281-298 (May 1961). 
[28] Young. H. and A. Kooharian, "A Decision Model for the Operational Acceptance Test (OAT) of Minuteman Ground Elec- 
tronic Equipment," Sylvania Electronic Systems, Res. Rept. No. 477 (Sept. 1965). 
[29] Zelen. M. (Ed.). Statistical Theory of Reliability, University of Wisconsin Press, Madison, 1963. 



GOAL PROGRAMMING IN ECONOMETRICS 

W. Allen Spivey 

University of Michigan 

and 

Hirokuni Tamura 

University of Washington 



1. INTRODUCTION 

Optimization models have been suggested as an alternative to simultaneous equation models for 
use in dealing with the problems of policy formulation and decision making at the level of the national 
economy (Theil [11]; Holt [5]; van Eijk and Sandee [14]; and Hickman [4]). In simultaneous equation 
models a policy maker prescribes a set of fixed values for goal or target variables and these are inter- 
related with a set of instrument variables by means of a system of simultaneous linear equations. The 
policy maker then solves this system (in which the number of equations or restraints is usually required 
to be equal to the number of unknowns) for the values of the instrument variables. It has been empha- 
sized (Tinbergen [12], [13]) that the policy maker's welfare function is not explicit in this model but is 
implicit in the selection of the goal and instrument variables. In the optimization approach one develops 
a welfare (preference) function which is optimized with respect to a set of econometric restraints. As 
is usual with such models, there can be many feasible solutions, and the optimization process then 
chooses a best or optimal solution from the set of feasible solutions. 

In the Theil optimization model (Theil [11]) a preference function is developed in which both goal 
and instrument variables appear. A desirable configuration of goal and instrument variables, not neces- 
sarily satisfying the constraints, is specified by the policy maker, and a quadratic function in deviations 
between the desired and feasible values of the variables is optimized (minimized). Although the quad- 
ratic function allows the policy maker to take account of a decreasing marginal rate of substitution 
between any two variables (at least over an appropriate subset of the domain of the function), it has the 
disadvantage of assigning the same penalty to overattainment of a goal as to underattainment (Hick- 
man [4], p. 6). A policy maker, for example, might not want to accept underattainment of a goal, but 
might be indifferent as to whether the goal was met or exceeded.! 

In this paper, we introduce an optimization model which differs from that of Theil in several ways. 
It permits overattainment and underattainment of any goal, yf, to be weighted either equally or dif- 
ferently in the decision maker's preference function. The model can be regarded as a special type of 
linear optimization model (as we make clear below) so it is easy to handle computationally and a wide 
range of parametric or sensitivity analyses can be implemented on a computer. These analyses can 



tStill another approach has been developed (Spivey and Tamura [9]), that of using a generalized simultaneous econometric 
model (GSEM), where the word generalized is used to emphasize that the number of equations or constraints is not required to 
be equal to the number of variables. The GSEM can be regarded as a link between conventional simultaneous equation models 
and the Theil model since, like the former it consists of simultaneous linear equations and like the latter a GSEM can be shown 
to have a solution which is identical with an optimal solution to a Theil type optimization model. 

183 



184 W. A. SPIVEY AND H. TAMURA 

display interactions between alternate optimal policies and goal attainment and can assist the policy 
maker in studying explicitly many of the trade offs between underattainment, attainment, and over- 
attainment of selected goals. Such sensitivity analyses are not available for either the Theil model 
or Tinbergen type simultaneous equation models. The model can also easily accommodate explicit 
restrictions on instrument variables in the form of nonnegativity and upper and lower bounds. 

The model is that of goal programming, originally developed by Charnes and Cooper in another 
context (see Charnes and Cooper [1]; also Charnes, Cooper, and Ijiri [2]; Ijiri [6]). We develop an appli- 
cation of this model using Klein's Model I of the United States economy as a basis (Klein [7]) and display 
in numerical form and interpret a sequence of decision making implications of several classes of the 
sensitivity analyses that can be performed with this model. 

2. GOAL PROGRAMMING AND GOAL ATTAINMENT 

A reduced form of a linear econometric model can be written 

(1) y=Rx + s, 

where y is an n by 1 vector of goal or target variables, x is an m by 1 vector of instruments,/? is ann by m 
(real) matrix of impact multipliers, and 5 is an n by 1 vector of constant or additive terms. If m — n, this 
is the reduced form of a conventional system of simultaneous linear equations; if m ^ n, then (1) is an 
example of a generalized simultaneous econometric model, in which m = n is a special case (Spivey and 
Tamura [9]). 

Let y* be the goal vector prescribed by the policy maker and consider the model 

(2) y* = /?* + s, 

where R and 5 are assumed to be known. If R is singular, then (2) either has no solution or infinitely 
many solutions. If R is non-singular, then (2) always has a unique solution; however, it may not be 
acceptable to the policy maker because one or more components violate restrictions which are present 
in the mind of the policy maker, but are not explicitly a part of (2) — for example, nonnegativity.t 
We reformulate (2) as 

(3) y*-s = Rx-Iz + + h- 

x,z + ,z~mO, 

where / denotes an n by n identity matrix, R is the n by m matrix introduced in (2), and z + and z~ are 
n by 1 (unknown) deviation or discrepancy vectors. The ith components zf of z + and zf of z" measure 
the deviation upwards and downwards, respectively, of feasible y, values from the corresponding goal 
values yf; in other words, the ith constraint of (3) can be written 

(4) yt-s^RfX-zt + zi (»=1, ...,»), 

where /?, denotes the ith row vector of the matrix/? and s, the ith component of the vectors. Note that 
a nonnegativity restraint on x is explicit in this formulation. 

If the policy maker can find a feasible instrument vector x for which bothz+ andz,~ in (4) are zero, 
then the goal yf can be attained exactly. On the other hand, if there is no feasible vector x for which 
zt—zf = in (4) he cannot attain the prespecified goal yf. He can, however, find a feasible x that will 



tlf implicit restrictions on x are violated in the case R is nonsingular, Tinbergen [12], [13] would suggest model reformulation 
by selecting another goal vector y*' in an ad hoc manner, substituting y*' for y* in (2). and again solving for x. 



GOAL PROGRAMMING 185 

allow him to come as close as possible to y* in a sense we will soon make clear. Now consider the 
following goal programming model 

(5) minimize e T z + + e T z~, 
subject to 

(6) y*-s = Rx-h + + Iz- 

x,z + ,z-^0, 

where e T is an n by 1 row vector t each of whose components is 1. 

The z+ and zf (i= 1, . . . , n) enter into the objective function (5) with a positive coefficient so the 
minimization process reduces the values of these variables in an optimal solution to the lowest feasible 
level and (5) attains a minimum value of zero if and only if both z + and z~ are the zero vectors (if and 
only if zt=Zi~ = Q for each i). The equations (6) clearly have feasible solutions because the rank of 
the matrix / assures that (6) can be solved for the z's alone. 

The advantage of reformulating (2) as a goal programming problem should now be clear: if (5) 
has a minimum value of zero, then (2) has a solution vector x which enables the policy maker to attain 
his prespecified goals yf (1 = 1, . . ., n) exactly. If it does not, then (2) is not solvable for a feasible 
instrument vector x. In the latter case, however, an optimal solution to (5) and (6) will disclose (a) 
how close the policy maker can get to the desired goal vector y*, and (b) whether this approximation 
is on the high side or low side for each individual goal, y* . 

The statement (b) can be made, and nonsense values of zf and zj avoided, because the vectors 
in the system of restraints (6) associated with zf and zj are negatives of each other. It is impossible 
for both zf and z,~ to assume positive values in an optimal solution because the basis property of 
linear independence is preserved in the simplex solution procedure. Denote an optimal solution to 
(5) and (6) as 

'Xo 

(7) 




where xj= (*i,o v • ., x m ,o), Zq T = (zt, o> • • •'> z n, o) an( ^ z o T= ( z i~, o» • • •> z n , o)- There are exactly three 
possibilities in an optimal solution to (5) and (6): 

(i) zt,o = Zi,o — 0, 

(ii) <•> 0,^0=0, 

( iH ) z<to=0, zf, >0. 



t The superscript T denotes a transpose. All vectors will be regarded as column vectors unless this superscript is used. The 
model (5) and (6) can also be written 

Minimize [0, e T , e r ] z* 
Subject to V 

y*-s= [/?,-/, /] lz< 

x 

z* \ 



186 W. A. SPIVEY AND H. TAMURA 

Case (i) means that the ith goal y* is attained exactly whereas cases (ii) and (hi) indicate that a deviation, 
upwards and downwards respectively, must be accepted by the policy maker desiring a nonnegative 
solution to (2). 

The loose phrase, "how close the policy maker can get to the goal vector y*," is given a precise 
meaning in goal programming because in (5) and (6) one is minimizing 

W X \yf-(Rix+ Si )\, 

i 

where the vertical bars denote absolute value. An optimal solution x yields values of the discrepancy 
variables Z+ and zj which have the following effects: 

(i') Zi,o = z7,o =Q if y* = RiXo + sr, 

(ii') zt >0 =RiXo+Si — yf, z f " =0 if y?<R(Xo+Si; 

(iii') Zi,o=y*— (Rtxo + Si), zl o =0 i( y*> RiXo + si. 

This can be stated in still another way. Given an xo in an optimal solution (7), the corresponding target 
configuration y , which occurs when xo is used in (2), is 

(9) y = Rx + s. 

If case (i) occurs, yo = y* as desired, but in either case (ii) or (iii) y ¥= y*. In the latter cases yo is the 
closest vector to y* where closeness is measured by the absolute value metric or distance measure 
(this metric is sometimes called the L t metric). That is, yo is an optimal solution to the problem 

n 

(10) minimize^ \y*~ Yi\, 

subject to 

(11) y = Rx + s,x^0. 

The measure of closeness can be generalized by using an L p metric and minimizing 

(12) (2ly<*-*| p V 

subject to the restraints (11), where p is any integer greater than or equal to one. If p = 2, this is the 
same as minimizing the quadratic function X(y* — y;) 2 subject to the restraints (11), which is the form 
of the objective function in the Theil optimization model. The model in this case selects a feasible 
vector y which is closest to y* in the sense of Euclidean distance or ordinary least squares. 

The goal programming model (5) and (6) can be augmented so as to contain political and institu- 
tional as well as nonnegativity restrictions on the instrument variables x,(i = l, . . ., m). Such 
restrictions, called boundary conditions by Tinbergen ([12], [13]), must be dealt with implicitly in 
the use of his models, but can be incorporated as upper and lower bounds on xi. These offer no difficulty 
in a linear programming model. For example, we could 

(5) minimize e T z + 4- e T z~ 

subject to 

(6') y* - 5 = R X - Iz + + Iz- 



GOAL PROGRAMMING 187 



Xi^ Li 


US, 


: k ^ U k 


keT, 


X, z + , z~ 


i=0, 



where S and T are any subsets of the set of integers {1, . . ., m} and Li and Uk are specified lower 
and upper bounds, respectively. The model (5) and (6') can be regarded as a decision model which 
explicitly takes account of all relevant restraints, both econometric and political or institutional. 

3. WEIGHTING DISCREPANCY VARIABLES 

In addition to determining a feasible vector y which is closest to the desired vector y* in terms 
of the L\ metric, a policy maker might regard attainment of certain goals as more important than that 
of others and might wish to use different weighting factors on deviation variables to reflect this. One 
can accordingly replace e T in the objective function (3) by other nonnegative weighting vectors u T 
and v T . This gives the following objective function instead of (5), 

(13) u T z + + v T z~ 
which is minimized subject to (6). Also, 

u T = v T = e T 

is a special case in which each deviation variable is given a weight of one. Thus if a policy maker 
desires to attain a given goal yf more than others he can attach larger weights m and v ,• to the corre- 
sponding deviations zf and z,~, respectively. If he would accept an overattainment of yf, but cannot 
tolerate any underattainment, he could set Uj = 0, but would make Vj large enough to prevent zj from 
assuming a positive value in an optimal solution. t In this case, the model will allow the policy maker 
to impose the following condition 

(14) yf =i R jX + s h 

i.e., we can accommodate an inequality restraint if desired — a feature which is not incorporated in 
the Theil optimization model. 

Moreover, one can parameterize the objective function in a goal programming model to see how 
an initial optimal deviation package — say, one in which some goals are attained exactly and some are 
not attained initially — varies as different nonnegative weights are attached to the deviation variables. 
This can be contrasted with the effect of a weighting factor in a Theil model (L 2 metric model), which 
permits a delicate control of the size of the corresponding deviation symmetrically above or below 
the goal vector y* (for an empirical study of the difference in effects of the use of Li metric and L 2 
metric, see Tamura [10]). 

4. SOME EXAMPLES OF DECISIONS IN GOAL PROGRAMMING MODELS 

The classical Klein Model I (Klein [7]) will be used to present some numerical examples of goal 
programming decisions. Model I in reduced form is shown below; actual values have been given to the 
lagged variables and the trend variable so that model can be used in examining 1933 policy decisions. tt 



tFor any nonnegative vector of coefficients [u T v T ], (13) has a lower bound of zero in a minimization problem. Only models 
with nonnegative coefficients in the objective function will be considered in this paper. 

ttThis model was also used for illustrative purposes by Theil [11] and Spivey and Tamura [9]. 



188 
Consider 

(15) 



W. A. SPJVEY AND H. TAMURA 




0.6659 


-0.1883 


0.6705' 


0.2238 


-1.2808 


1.1188 


0.1620 


-0.2033 


0.8102 


0.0517 


-0.2959 


0.2584 


0.6142 


-1.4842 


1.9289 


0.0517 


-0.2959 


0.2584 




+ 



/37.58^ 
' 9.0764 

23.0109 

-5.4938 

32.0873 
i 201.6060/ 



where C= consumption expenditures, P= corporate profits, W\ — private wages, /= private invest- 
ment, Y— national income, and K= end of year capital stock are endogeneous variables and W<l = gov- 
ernment wages, 7"= tax revenue, and G = government expenditures are exogeneous variables (all 
variables are measured in terms of billions of 1934 dollars). The constant terms are calculated under 
the assumption that the random disturbances are forecast accurately. 

EXAMPLE 1. Parameterization of Weights: Suppose that the goal vector selected by a policy 
maker is y T = (C, P, Wi, I) and the instrument vector is x T = {W2, T, G,). Then in the manner of Tin- 
bergen and Theil we can formulate a linear econometric model in these vectors y and x by extracting 
the first four equations from the reduced form (15). 

Let the goal values be set as a uniform 10 percent increase over the actual values of 1932 for each 
goal variable, i.e., y* T — (50.16, 7.7, 31.9, —5.58). Then we use the following goal programming model 
to examine 1933 policy decisions, 



(16) 
subject to 

(17) 



minimize zt + z c + zp~ + z p + z^ x + z Wl + zf + Zj 




+ 



0.6659X 
-0.2238 \ w 
-0.1620 I 
-0.0517/ 

-i\ , / 1 





0/ 




0.6705\ 
1.1188 G + 
0.8102 



0.2584 






U*i + 




W 2 , l , tr, Z c , Z c , Z p , Z p , Z Wl , Z Wi , Z/ , 2/ = U, 

where (16) and (17) correspond to (5) and (6), and where each goal is given equal weight. An optimal 
solution to this problem leads to the policy 



(18) 



6.870 1 
I 13.650|, 
,15.770/ 



and the corresponding target configuration yo is 



(19) 



yo- 




GOAL PROGRAMMING 



189 



Comparison of yo with the prescribed vector y* indicates that the three goals C*, P*, and W* are 
attained exactly and the private investment goal /* is underattained by $.23 billion. 

Suppose that the policy maker, upon being informed of (19), wants to assign a larger weight to the 
attainment of the investment goal/* than is the case in the objective (16). Then the deviation variable 
zj can be given the weight A (z/ being given the same weight 1), and we 



(20) 



minimize z£ + z c + zp~ + z p + z+ 1 + z Wi + z/" + kz, 



subject to (17). We study this optimization for various values of k (continuous parameterization of 
k). Table 1 summarizes relevant information. For 1 Si k < 4.328 an optimal solution remains unchanged 
and for each such choice of weighting factor k, Xo in (18) is optimal and y in (19) is realized. This is indi- 
cated by the first data row of Table 1. For weight k ^ 4.328, we must shift to a new optimal solution x , 
getting a corresponding y . The second data row in Table 1 displays values of these vectors. The invest- 
ment goal is attained for weights k in the range indicated and we overattain the goal P* by one billion. 

Table 1. Optimal Policies and Consequent Goal Attainments 



Range of Values 
of A 


Components of y 


Components of xo 


Co 


Po 


r,, 


h 


w 2 , 


T 


Co 


lgK 4.328 


50.16 


LI 


31.9 


-5.71 


6.870 


13.650 


15.770 


4.328 S k < °° 


50.16 


8.7 


31.9 


-5.48 


6.845 


12.641 


15.512 



Note: Underlining indicates exact attainment of goal. 

Comparison of the two optimal policies indicates that attainment of /* is brought about by in- 
creased deficit spending; Go - To = 15.770— 13.650 = 2.120 billion for the original optimal x , whereas 
G — To = 15. 512— 12.641 = 2.871 for the alternate optimal xo- Note that the increased deficit is not 
attributable to increased values of Go, but rather to decreased values of T . That is, in this model 
optimal adjustments call for reducing taxes rather than increasing government expenditures. 

A policy maker might wish to inform himself further about the interactions between T, G, I, and 
optimal policies. T, G, and / are related by means of the fourth equation in (15), 



/ = - 0.0517r 2 - 0.2959r+ 0.2584G - 5.4938. 

One could study how sensitive an optimal policy is to changes in either of the multipliers — 0.2959 and 
0.2584. A computer program can be used to parameterize the matrix in (17), one row or one column at 
a time, so parameterizing any one element in the matrix is a special case. Specifically, one could 
determine a range of values containing — 0.2959 within which an xo remains optimal, assuming that 
other columns of the matrix are held constant. Similar remarks can be made concerning the para- 
meterization of the government expenditures multiplier 0.2584. Indeed, the following kinds of para- 
metric analysis can be performed on the model (16) and (17): (a) nonnegative parameterizing of the 
coefficients in (16); (b) parameterizing the column vector of constants in (17); (c) parameterizations of 
(a) and (b) simultaneously; (d) parameterizing the matrix in (17) one row or one column at a time. Thus 
the model provides the policy maker with considerable flexibility in both policy formulation and decision 
making. 



190 W. A. SPIVEY AND H. TAMURA 

EXAMPLE 2: Prespecification of Both Goals and Policies: Theil [11] introduced an optimization 
model in which the policy maker specifies a goal vector, y*, and a desired instrument vector, x*. We 
reformulate this model by means of goal programming and introduce additional flexibilities by using 
alternate objective functions. 

The target or goal variables chosen by the policy maker are C, I, and D, where D is a distribution 
variable defined as W x — 2P (see Theil [11], p. 79). Prespecified goal values are, in billions of 1933 
dollars, C* = 49.69, /* = — 3.10, and D* = 11.25. The desired instrument values for W 2 , T, and G are 
5.04, 7.40, and 10.44, respectively (for specification of goal and instrument variables, see Theil [11]). 
Thus 

(21) y* T = (C*,I*,D*) = (49.69 -3.10 11.25) 
and 

(22) x* T = {Wt,T*,G*) = (5.04 7.40 10.44). 
The goal programming model is to 

(23) minimize e T z + + e T z~ 
subject to 

(24) /y*\ /s\ (R 




X, Z + , 2^0, 

, (A . 

where I in (24) comes from (21) and (22), and where 

/ 37.5812 \ / 0.6659 -0.1833 0.6705' 

s = -5.4939 I, R= -0.0517 -0.2959 0.2584 

\ 4.8580/ \ 0.2855 2.3583 - 1.4274, 

were calculated from data in Klein's model 1 in the manner suggested by Theil (see Theil [11] p. 87 
and Tamura [10], p. 46), h and h denote identity matrices of size 3 and 6, respectively, and 

Z + = (Zf , Zj , z D , z W2 , Z T , Z (: ) , 
z~ = {zj}, ZJ, zjj, z W2 , Z T , z a ) , 
x T = (W 2 , T,G). 

As observed earlier, the objective function we are using assigns the same weights to each zj 1 " and 
zf. If the policy maker would accept an overattainment of goals C* and /* and an underattainment of 
goal D*, then minimization of the following objective function subject to (24) could be suggested, 



(25) f 



where/ is a column vector which contains 12 elements all but 3 of which are set equal to 1. The remain- 
ing three elements are the coefficients of zf , zf, and zj, and are set equal to zero. The zero coefficient 
of these variables then essentially excludes them from the minimization process. 



GOAL PROGRAMMING 



191 



Optimal policy solutions under the two different objectives are given in Table 2. For purposes of 
comparison we have also indicated a Theil type optimization decision, i.e., one from a L 2 metric model. 

TABLE 2. Alternate Optimal Policy Decisions (Billions of 1934 Dollars) 





Values of y 


Values of x 


C 


I 


D 


w 2 


T 


G 


Desired Values 


49.69 


-3.10 


11.25 


5.04 


7.40 


10.44 


Optimal Solution 
with (23) 


46.35 


-5.55 


11.25 


5.04 


8.42 


10.44 


Optimal Solution 
with (25) 


46.54 


-5.25 


8.84 


5.04 


7.40 


10.44 


Theil Model 


47.77 


-5.41 


10.42 


6.44 


8.32 


11.13 



Note: Underlining indicates exact attainment of goal. 

5. CONCLUDING REMARKS 

Earlier we indicated a generalization of optimization models in terms of the L p metric, a generali- 
zation which includes both goal programming and Theil type optimization models as special cases, i.e., 



subject to 



minimize I ^ \yf — J\ | p ) 1/p , 
y = Rx + s, x ^ 0. 



Although for p = 2 this is a nonlinear optimization problem and, as such, may appear tQ present 
insurmountable computational difficulties, recent research (Duffin, Peterson, and Zener [3] and 
Peterson and Ecker [8] ) has shown that this problem is one of a class of nonlinear optimization prob- 
lems known as geometric programming. Results for these models are available whereby numerical 
calculations and approximations can be developed. Moreover, certain types of duality theorems and 
"complementary slackness" properties have been established, the economic interpretations of which 
are yet to be explored. Research by the authors on these matters is now in progress. 

REFERENCES 



[1] Charnes, A. and W. W. Cooper, Management Models and Industrial Application of Linear Programming (John Wiley & 

Sons, Inc., New York, 1961), Vol. 1. 
[2] Charnes, A., W. W. Cooper, and Y. Ijiri, "Breakeven Budgeting and Programming to Goals," J. Accounting Research /, 10-43 

(1963). 
[3] Duffin, R. J., E. L. Peterson, and C. Zener, Geometric Programming— Theory and Application (John Wiley & Sons, Inc., 

New York, 1967). 
[4] Hickman, B. G. (editor). Quantitative Planning of Economic Policy (The Brookings Institute, Washington, D.C., 1965), 

Chap. 1. 
[5] Holt, C. C, "Quantitative Decision Analysis and National Policy: How Can We Bridge the Gap?" in Quantitative Planning 

of Economic Policy, edited by B. G. Hickman (The Brookings Institute, Washington, D.C.. 1965). 
[6] Ijiri, Y., Management Goals and 'Accounting for Control (North Holland Publishing Company, Amsterdam; Rand McNally 

Company, Chicago, 1965). 
[7] Klein, L. R., Economic Fluctuations in the United States 1921-1941, Cowles Commission Monograph No. 11 (John Wiley & 

Sons, Inc., New York; Chapman and Hall, Ltd., London, 1950). 



192 W. A. SPIVEY AND H. TAMURA 

[8] Peterson, E. L., and J. G. Ecker, Geometric Programming: Duality in Quadratic Programming and L v — Approximation, 

The University of Michigan, Technical Report No. 1968-2, 1968. 
[9] Spivey, W. A., and H. Tamura, "Generalized Simultaneous Equation Models." to be published in International Economic 

Review. 
[10] Tamura, H.. "Linear Models for Macroeconomic Policy Making" (unpublished Ph.D. thesis) The University of Michigan, 

1967. 
[11] Theil, H., Optimal Decision Rules for Government and Industry (Rand McNally and Company, Chicago, 1964). 
[12] Tinbergen, J., "On the Theory of Economic Policy," of Contributions to Economic Analysis (North Holland Publishing 

Company, Amsterdam, 1952), 2nd ed., Vol. I. 
[13] Tinbergen, J., Economic Policy: Principles and Design of Contributions to Economic Analysis (North Holland Publishing 

Company, Amsterdam), Vol. XI. 
[14] van Eijk, C. J. and J. Sandee, "Quantitative Determination of an Optimum Economic Policy." Econometrica 27, 1-13 (1957). 



PRODUCTION AND EMPLOYMENT SCHEDULING IN MULTISTAGE 

PRODUCTION SYSTEMS 

Christoph Haehling von Lanzenauer* 

The University of Western Ontario 
London, Canada 

ABSTRACT 

A model is developed for planning optimal production and employment levels in multi- 
product, multistage production systems. The market requirements for each product over 
the planning period are given. Backorders and/or shortages are permitted. Backorders 
and shortages must be considered in order to determine the amount of each product's demand 
that should be filled, backlogged, or go unsatisfied if the production capacity is insufficient 
to fill all market requirements. Backorders and shortages, on the other hand, are desirable 
under certain dynamic market conditions. 

I. INTRODUCTION 

The production and employment scheduling problem — planning production and employment levels 
over time — has received a great deal of attention. Models and decision-rules generally based on single- 
stage production systems have been developed for many situations under various conditions [7]. When- 
ever many products are to be manufactured the concept of aggregate production scheduling has 
frequently been suggested [1, 5]. The aggregate production level is then represented by some appro- 
priate measure of production per unit of time [5] in view of the common denominator required. Once 
the aggregate production level is specified the problem of determining the quantities and the timing 
for the individual products still remains. This problem of disaggregation is of special importance when- 
ever the available capacity is insufficient to produce all market requirements and therefore some 
demand has to go unsatisfied. It is not at all obvious how the disaggregation decisions should be made 
in view of insufficient capacity. 

Following along the lines of more recent work in production planning [2, 3, 8], we consider a multi- 
stage production system consisting of various production facilities (stages) linked together by in-process 
inventories. The objective of this paper is to develop a model for multiproduct, multistage production 
systems by which optimal production decisions on a disaggregate basis for both sufficient and insufficient 
production capacity can be made. As a point of departure we take Zangwill's [8] model. The subsequent 
analysis, however, differs from Zangwill's work on a number of points: 

(a) The demand requirements need not be met. 

(b) Each production stage manufactures various different products. 

(c) Workforce decisions are integrated. 

(d) A different cost structure is assumed. 

II. THE MODEL 

The Production System 

The elements of a multistage production system are the individual production stages consisting 
of production facilities and an. inventory. These production facilities can be thought of as multipurpose 



* Assistant Professor of Business, School of Business Administration. 

193 



194 



C. H. von LANZENAUER 



machinery. The output of a given production stage will be fed into the inventory bank from which the 
units can be used either as inputs to higher numbered production stages or to supply market require- 
ments for the stage's products. Figure 1 describes such a production system. 



"ff~ 



MARKET 



Production and Inventory Aspects 

Let Dj, be the market requirements of production stage fs (j= 1, 2, . . . , J) product i (i= 1, 2, 
...,/) during the period t(t=l, 2, . . ., T). It is assumed that the requirements in each period 
over the planning period are known in advance. The planning period is medium range so no discounting 
is required. The market requirements are expressed in units. Note that the index i at each production 
stage represents physically different products. Let Hj, be the amount of product i in the inventory of 
stagey at the end of period t. The inventory equation can then be expressed by (1). 



(1) 



tih = H h 



+ X ( +X j 



D h 



k=j+l n=l 



for all i, t andy= 1, 



,y-i 



with XjtR {Xj l0 ) being the amount of product i produced in stage j during period t on regular time 
(overtime). aj£ (ajg 3= 0) can be considered as a production coefficient and is the number of units of 
production stager's product i required to produce one unit of stage &'s (k=j+l, j+2, . . ., J) 
product n (n= 1, 2, . . ., /). aj£ = O if product i of stagey is not being used for the production of 
stage Fs product n. The inventory equation for the last stage J of the production system is 



(2) 



Ui — JJi _l_ Vi _L Yi _ f)i 



for all i, t. 



The formulation of (1) and (2) assumes that the production capacity available is sufficient to satisfy 
all market requirements. If market requirements exceed the available production capacity, the excess 
demand must either be backlogged into the next period or has to go unsatisfied. The possibility for 
backlogs and shortages is only meaningful for market requirements since the requirements of sub- 
sequent production stages, i.e., the production levels, are variable themselves. The inventory equation, 
adjusted by backlogs and shortages, can then be expressed by (3) 

(3) n}, = Hj t _ l +X} at +X} lo -(Dl t + B} t _ l )+Bj t +Sl t - ^ z, -» *"*m > "tto 



k=j+] n=i 



for all i , t and y = (/' + 1 , 



J -I) 



with #j ( being the demand backlogged into period £+ 1 and Sj, the amount of the periods total demand 
(Dj t + Z?j,_,) which will never be filled. Note that the variables//],, Bj t , and SL are nonnegative vari- 
ables. Generally, one can assume that, whenever a positive inventory exists of production stage y's 
product i at the end of period t , Bj, and Sj, will be zero; on the other hand, whenever #j, and/or Sj, is 
greater than zero, there will normally be no inventory available. It should be realized, however, that 
the formulation of (3) permits the possibility of having positive inventory as well as backlogged and/or 
unfilled demand at the same time. This configuration can be desirable considering over time changing 



MULTISTAGE PRODUCTION SYSTEMS 195 

economic conditions. For example, in view of a speculative increase in the selling price in periods 
ahead, it might be advantageous to let some of the current periods demands go unfilled, though inventory 
is available. The inventory equation for the last production stage is 

(4) //', = //),_, + X> m + X* M - (/>*,+ 0j,_, ) +B< ( + Sj, , for all i, t. 

For the special case of a production line (i.e., production stager's output will only be used in produc- 
tion stage;' + 1) we have instead of (3), 

(3a) ^ = ^- 1 +^-^o-(^ + ^- 1 )+^ + ^-2 «&, C*£„»+*W' 

n=l 

for aU i, t. 
In view of hmited in-process inventory capacities, we can restrict 

(5) EP/ffJ* 1 * for all;,*, 

i=\ 

where /3j is the amount of space required to store one unit of production stager's product i. Lj t repre- 
sents the available inventory capacity in space terms. Correspondingly, restrictions can be placed on 
the maximum amount to be backlogged per period. Such an upper limit is desirable for Bj T in view 
of the fact that the cost associated with backlogs are normally less than those associated with shortages. 
The production quantities are subject to the following capacity restrictions 

(6) X'W** 5 ^ for all;,, 
and 

(7) S 7 /*^ * for all;,,, 

where t) is the time required to produce one unit of production stage ;'s product i and Rjt (Ojt), the 
capacity in time units of regular time (overtime) available on stage; during period t. 

Workforce Considerations 

The production decisions can only be carried out if the necessary workforce is available. The 
workforce decisions are to be made in view of the production requirements and with respect to the 
dynamic costs involved. (For a discussion of dynamic costs see [5]). The workforce available during 
period t will be denoted by W t . According to [6] we can differentiate between "productive" (Wf) 
and "idle" workforce {W}). So 

This distinction is desirable in view of the possibility of keeping workers on a reduced pay rate rathei 
than laying them off. Let y) be the number of units of stage ;'s product i which can be produced bj 
one worker during period t. Assuming a constant worker productivity (i.e., a constant yj), the require< 
productive workforce for regular time or overtime production in t can be expressed by 

(9) ^-is^ii 



196 C. H. von LANZENAUER 

and 

<io) ri-22$*ir 

The productive workforce in period t is therefore 

(id pf=£i^c^+*jw- 

The total workforce available in period £ can then be expressed by 

(12) £i^(^+^) + ^*i£^(^-«+^-io) + »7-i + ft-»-^-i for all*, 

with G?_i being the number of workers hired at the end of period t— 1 and Ff_i the number of workers 
laid off at the end of t— 1. Since the decision of which workers to lay off is usually made on grounds of 
seniority, we don't have to specifically account for this fact within the model. The same holds for the 
"idle" workforce which is generally determined by juniority. The number of workers to be assigned to 
the various production stages during period t can easily be obtained from (9) and (10). 

From an organizational point of view, it is necessary to maintain the proper relation between the 
number of processing workers and supervisory labor. Let Wf be the number of supervisory labor in 
period t, then we have an analogous situation to (12) 

(13) W t s =W t s . 1 + M l .i-N t - 1 for all f, 

with M t -\ being the number of supervisory labor added and N t -\ the number laid off at the end of 
period t— 1. The relation between the processing and supervisory labor can then be expressed as 

(i4) \t±j j V!m+xj»)'*rf. 



j=l i- 



with 8 being the span of control. 



The Objective Function 

The objective of the model is to determine the production and employment levels so that the con- 
tribution to profits over the planning period is maximized. 

Let r'j t be the selling price of one unit of production stager's product i during t. If we assume that r' 
is constant for all t , the total revenues of the production system are 

iiirjiDj-sj,), 
t = i j = i i=i 

since any amount backlogged into a future period will bring the same revenue. Assuming an over time 
changing selling price rj t , the total revenue depends on the actual price any backlogged demand can 
ultimately be sold for. If backlogged demand Bj t has to be filled at the existing price rk of period t the 
total revenues are 

t=i j=i i=i 



MULTISTAGE PRODUCTION SYSTEMS 197 

If on the other hand backlogged demand will be filled at current prices, the total revenues become 

t=i }=i i=i 

Production costs are the sum of material and labor costs. Capital cost will be considered as fixed. The 
cost to produce stager's product i depends not only on the decisions made on stagey but also on the 
production and employment decisions made on the lower-numbered stages which provide input to 
stagey. However, if one handles the labor costs en bloc, the material costs of stagey 's product i are only 
the costs for additional raw material since the costs of the inputs from lower numbered stages are 
already accounted for. Let cj t be the cost of additional raw material to produce one unit of stage y's 
product i during period t. The total material costs are then 

t=\ j=i 1=1 

The inventory costs are the product of the average amount in inventory and an inventory cost factor 

z). Therefore, we have 

j 

t=i j=i i=i ^ 

Backlog and shortage costs are 

t=\ j=i i=i 

with bf being the penalty for backlogging one demand unit of stage /s product i for one period, s! is 

is the expected cost in terms of loss of goodwill, over and above the lost revenue, of one unit of stage 
/s product i that cannot be filled. 

The costs associated with the labor-force are 

t =i j=i ,=i J 

with p being the regular-time pay rate of a processing worker per period, k the overtime premium, 
and A the pay reduction for the idle workforce. h(hs) is the hiring cost per worker (supervisor), f(fs) 
the cost of laying off one worker (supervisor), and ps the supervisor's pay rate per period. Combining 
the above expressions* leads to the objective function (15) 



/ 

I 

t=\ j=\ i=\ 



(is) tp=z j 2 W-sp-^«+*y 






"Assuming an overtime constant price r' . 



198 C. H. von LANZENAUER 

- hG, ~/F t - p s Wf - h s Gf -fsFf} 

-iii{«-(4+^-(4+*^. 

z=i j=\ i=i "■ N ry \ Yj ' 

-- z\{H\ +H i l )-b i .B i - (r?+sf)S* 

-p{\-\)Wl-h-Gt-f-F t -p s W s t -h s Mt-f s Nt}- 

Maximizing (15) with respect to (3) or (3a), (4), (5), (6), (7), (12), (13), (14), and the nonnegativity constraints 
for the variables is a problem that can be solved by linear programming.t 

III. CONCLUSION 

The purpose of model was to formulate the production and employment scheduling problem 
on a disaggregate basis in multiproduct, multistage production systems. The model determines the 
amount of the demand for each product that should be satisfied, be backlogged, or remain unfilled. 
These results are important whenever the production capacity is insufficient to produce all market 
requirements. On the other hand, the formulation can be very useful under certain dynamic market 
conditions. The model determines also how the productive workforce is to be assigned to the various 
production stages. The model seems to be applicable to integrated iron and steel manufacturing systems 
as well as to various chemical processes. 

The author wishes to thank Andrew Grindlay and Howard E. Thompson for valuable comments and 
suggestions. 

REFERENCES 

f 11 Buffa, E. S., Operations Management — Problems and Models (John Wiley & Sons, Inc., New York, 1968) 

[2] Elmaghraby S. E. and A. S. Ginsberg, "A Dynamic Model for the Optimal Loading of Linear Multi-Operations Shops," 

Management Technology, Vol. 4, No. 1 (1964). 
[3] Haehling von Lanzenauer, C, "Formulation of the Production Scheduling and Sequencing Problem by Bivalent Linear 

Programming," Management Science, forthcoming. 
[4] Hanssmann, F. and S. W. Hess, A Linear Programming Approach to Production and Employment Scheduling, Management 

Technology, Vol., No. 1 (1960). 
[5] Holt, C, F. Modigliani, J. Muth, and H. A. Simon, Planning Production, Inventories, and Work Force (Prentice-Hall, Englewood 

Cliffs, N.J., 1960). 
[6] Lippman, S. A., A. J. Rolfe, H. M. Wagner, and J. S. C. Yuan, "Optimal Production Scheduling and Employment Smoothing 

with Deterministic Demands," Management Science, Vol. 14, No. 3 (1967). 
[7] Silver, E. A., "A Tutorial on Production Smoothing and Workface Balancing," Operations Research, Vol. 15, No. 6 (1967). 
[8] Zangwill, W. I., "A Deterministic Multiproduct-Multifacility Production and Inventory Model," Operations Research, Vol. 14, 

No. 3 (1966). 



tNote that ^ X 2 ri P\t is a con stant. 



A MUTUAL PRIMAL-DUAL LINEAR PROGRAMMING ALGORITHM 

Milton Y. Harris 
Naval Research Laboratory 

ABSTRACT 

A new primal-dual linear programming algorithm is exhibited. A proof is given that 
optimal solutions to both primal and dual problems (when such solutions exist) are found in 
a finite number of steps by this algorithm. A numerical example is included to illustrate 
the method. 

1.0 INTRODUCTION 

Since the publication of the primal-dual algorithm by Dantzig, Ford, and Fulkerson, in 1956 (see 
Refs. [3], [5], [7]), apparently little work has been done in this area (see Gomory and Balinski [1]). This 
paper presents a primal-dual algorithm which is similar to that of Dantzig, Ford, and Fulkerson in that 
it uses the complementary slackness concept and alternates between primal and dual problems; how- 
ever, it differs in that it is simpler, but requires a separate Phase I stage.* 

Essentially, this algorithm consists in generating and solving a sequence of subproblems which 
use the constraints of the original problem. The objective function at each stage, however, is determined 
by the solution of the problem of the previous stage. Whether this algorithm is any more efficient than, 
say, the revised simplex algorithm remains to be seen.t 

Before exhibiting the algorithm, we first establish some notational conventions and state some 
necessary results. 

2.0 NOTATION AND BASIC THEOREMS 

For the sake of simplicity, and without any loss of essential generality, we will consider only 

problems of the form 

maximize z—c'x 

subject to 

Ax^b 
(2.1) 

x^O 

where A is an {nxm) matrix, b is an {nx\) vector with nonnegative entries, c is an (mxl) vector of 
criterion function coefficients, and x is an {mxl) vector of variables. The "primes" attached to symbols 
denote the transpose of the indicated matrix. All entries in /I, b, and c are taken to be given constants. 
We may state the dual to this problem as 

(2.2) minimize u—w' b 



*For an explanation of the two stages of the simplex method, see Dantzig [4], ch. 5. 

t A Fortran program has been written for the CDC 3800 computer to compare the efficiency of this algorithm with the sim- 
plex method, but as yet, no conclusive results have been obtained. 

199 



200 M. Y. HARRIS 

subject to 

w'A S= c' 
(2.3) 

w' ^0 

where w is an (nxl) vector of variables. 

These two problems may be rewritten by introducing slack variables 5 and v for the primal and 
dual, respectively, and taking the transpose of 2.3 as follows: 

maximize z — c'x minimize u — w' b 

subject to subject to 

Ax + Is = b A'w — Iv = c 

(2.4) (2.5) 

{x,s)^0 {w,v)^0, 

where / is an appropriately dimensioned identity matrix in each case. 

We will digress here slightly to define formally the notion of complementary slackness, a concept 
which is central to the algorithm. 

DEFINITION: For given basic feasible solutions (x, s) to the primal and (w, v) to the dual, 
complementary slackness is satisfied if, and only if, 

i) w • s — 
and 

ii) x - v — 0. 

We conclude this section by stating four well-known results which are essential to our algorithm. 
Proofs are given by Dantzig in [4]. 

Suppose we have any two basic feasible solutions (x, s) for the primal (2.4) and (w, v) for the 
dual (2.5). Then 

LEMMA 1: w'b^c'x. 

LEMMA 2: (x, s) and (w, v) are optimal solutions for the primal and dual, respectively, if, and 
only if, w'b = c'x. 

THEOREM 1: (x, s) and (w, v) are optimal solutions for the primal and dual, respectively, if, and 
only if," complementary slackness is satisfied for these solutions. 

FUNDAMENTAL DUALITY THEOREM: If a feasible solution to 2.4 exists and z of 2.4 is bounded 
above, then optimal feasible solutions exist for both the primal (2.4) and the dual (2.5). 

3.0 THE ALGORITHM 

The algorithm is initiated with basic feasible solutions to both the primal and dual, either or both 
of which may be contained in the intital formulation of the problem or may be generated by Phase I of 
the simplex algorithm. We also assume that z of 2.4 has a finite upper bound. The fundamental duality 
theorem assures us that optimal solutions for both the primal and dual exist. 

Recall the primal and dual constraints may be written as 

(3.1) Ax + Is = b (primal) 
and 

(3.2) w'A-v'I=c' (dual). 



MUTUAL PRIMAL-DUAL ALGORITHM . 201 

Multiplying 3.1 by w' on the left and 3.2 by x on the right and applying Lemma 1 gives 

(3 3) ws + vx = wb — c-x^0. 

Let a = w • s + v • x. Now any decrease in the value of a brings us closer to the optimum since a decrease 
in a represents a movement toward the satisfaction of complementary slackness and a convergence of 
wb and c-x toward their common optimal value (see Lemma 2 and Theorem 1). 

We now proceed to define the algorithm. Suppose we have a basic feasible solution {x°, s°) to the 
primal and a basic feasible solution (w°, v°) to the dual with 

(x°, 5°) • {v°, w°)=x -v + s°-w > 

(i.e., the solutions are non-optimal). Let Y° be a current full-tableau matrix for the primal, F° be a cur- 
rent constant column for the primal, Z 1 be a current full-tableau matrix for the dual, and Z* be a current 
constant column for the dual. Consider the linear programming problem 

(3.4) minimize a — (x, s) • (v°, w°) 
subject to 

(3.5) Y°(x,s) = Y° 

(x, s) ^ 0. 
Suppose the optimal (if one exists) solution to 3.4 and 3.5 is (x 1 , s 1 ) and 

(zSs 1 )-^ ,^ )-^. 

Let Y 1 and YJ be the matrix and constant column, respectively, corresponding to this solution. If 
6*^ = 0, then {x 1 , s 1 ) and (w°, v°) are optimal solutions to the primal and dual, respectively, by Theorem 
1. If oto > (recall from 3.3 that a<f 3= 0), consider the problem 

(3.6) minimize a^ (x l , s 1 ) • (v, w) 
subject to 

(3.7) Z*{w i v)=Z\ 

(W, V) 52=0. 

Let the optimal tableau to this problem be represented by Z 2 and Z\ and the optimal solution be 
(w 1 , v 1 ); let 

(x^s 1 ) ■ (v 1 ,w 1 )=at 

As before, if af = 0, we are through. If not, we return to the primal and solve 

(3.8) minimize a 2 = (x, s) • (v 1 , w l ) 
subject to 

(3.9) ° 

(x,s) >0 

for (x 2 , s 2 ), Y 2 , Y§, and a* = (x 2 , s 2 ) ■ (v 1 , w x ). The algorithm consists in following this procedure of 
using the optimal solution at each stage as the criterion function coefficients for the next stage until 
for some n, a* = 0. 



202 M. Y. HARRIS 

In the next section we will show that the algorithm results in optimal solutions for both the primal 
and dual, whenever they exist, in a finite number of iterations. 

4.0 PROOF OF CONVERGENCE 

We may state the subproblems described above in general as follows: 

Primal 

minimize ak = (x, s) ■ (v^ k \ w [k] ) 

subject to 

YM(x,s) = YM 
(4.1) 

(*, s)2=0, 

where [k] — k/2 for A; = 0, 2, 4, . . . and 

Dual 

minimize a m = (x lm \ s [m] ) ■ {v, w) 

subject to 

ZM(w, v)=ZW 
(4.2) 

(w, v) 5=0, 

where [m] = (m— l)/2+ 1 for m — 1, 3, 5, . . . We first show that the subproblems 4.1 and 4.2 have 
optimal solutions. 

THEOREM 2: Both problems (4.1 and 4.2 above) have optimal solutions for any ke{0, 2, 4, . . .} 
and any me{l, 3, 5, . . .}. 

PROOF: Both problems (4.1 and 4.2) have feasible solutions by assumption. Now for every basic 
feasible solution (x, s) to the primal and (w, v) to the dual, we have by Lemma 1 

(4.3) x- v + w s^O. 

Since the problems in question are both minimizing problems whose objective functions are of the 
form of 4.3, we have just seen that both are bounded below by zero. By the fundamental duality theorem 
then, both have optimal solutions. 

Finally, we show that the algorithm does indeed converge. 

THEOREM 3: a„* = for some (finite) n. 

PROOF: We assume that the original problem has an optimal solution. Clearly 

a* < ao, 

since at worst (w°, v°) satisfies 3.7 and (x 1 , s 1 ) • (v°, w°) — c*o (this fact is just a general property of 
constrained optima, that is, that the optimal solution to a minimizing problem is no larger than any 
admissable solution). Suppose a* = a*. Then for every basic feasible solution (w, v) for the dual 

w • b — c • x x = w • s x + s l + x l • v (by 3.3) 
= «i? a* = ctQ (assumption) 

= w° • s 1 + x 1 ■ v° 
= w°-b-c-x 1 (by 3.3). 



MUTUAL PRIMAL-DUAL ALGORITHM 203 

Thus for every basic feasible solution (w, v) for the dual 

w • b — c • x 1 2= w° • b — c • x 1 or w • b 2= w° • b, 

which implies that (w°, v°) is the optimal solution for the dual. Similarly, for every basic feasible 
solution (x, s) for the primal 

w° • b — c • x = w° • s + x • v° (by 3.3) 

= «i 2= a* — ct * (assumption) 

= w°-b-c-x 1 (by 3.3) 

Thus for every basic feasible solution (x, s) for the primal 

w° • b — c • x ^ w° • b — c ■ x l or c • x ^ c ■ x 1 

which implies that (x l , s 1 ) is optimal for the primal. Now if a* > 0, then complementary slackness 
is not satisfied by the solutions (*', s 1 ) and (w°, v°) , and so at least one must be non-optimal by Theorem 
1. This contradiction (i.e., we just showed that if a* — a*, both (m; , v°) and (x 1 , s 1 ) are optimal) im- 
plies that if a* > 0, then 

a* < a *. 

We may repeat this argument for ai and oc-z and so on. We conclude therefore that 

a * > implies that a* +1 < a* for each A; = 0, 1,2,. . .. 

Since there are only a finite number of basic feasible solutions to the primal and dual, and since 
each iteration of the algorithm results in a strict decrease in a, if optimal solutions exist, then by 
Theorem 1, there exists an n (finite) such that 



Furthermore the solution to the problem for which a n is the objective function is the optimal solution to 
either the primal or dual with the solution to the preceeding problem being the optimal solution to the 
other. 

5.0 NUMERICAL EXAMPLE 

Consider the following primal problem (this problem was taken from Hadley [5], p. 269): 

maximize z= 2x x + 13* 2 + Tx 3 + 6x4 

subject to 

1 & 4*i + lx 2 + 3*3 + 5x 4 
3 3= — xi + 3x2 + 3x 3 — 3x 4 

(5.1) 7 3*2*1- 4x 2 +lx 3 + 8x4 

5 5s 3xi + 5x 2 — 2x 3 — 4x 4 

2 2= — 5xi — lx 2 + lx 3 + 6X4 
8 5=6xi + 2x2-9x3+lx4 



The dual is then 



minimize u = W\ + 3k;2 + 7w 3 + 5w 4 + 2w 5 + 8w 6 



204 



M. Y. HARRIS 



subject to 



(5.2) 



2 =£ 4wi — 1^2 + 2tt>3 + 3u>4 — 5w$ 4- 6wa 

13 =S 1 w\ + 3z#2 — 4m; 3 4- 5t^4 — \w$ + 2we 

7 =S 3wi + 7w 2 + Iwa — 2w4 + lws — 9w 6 

6 =£ 5^i — 3^3 + 8w 3 — 4w 4 + 6^5 + 1m;6 

After appending slack variables S\ through se to 5.1 and v\ through y 4 to 5.2 and writing both systems 
in detached coefficient form, we have 

Primal 



(5.3) 



b 


Xi 


*2 


x 3 


*4 


Sl 


«2 


53 


54 


55 


56 





2 


13 


7 


6 




















1 


4 


1 


3 


5 


1 

















3 


-1 


3 


7 


-3 





1 














7 


2 


-4 


1 


8 








1 











5 


3 


5 


-2 


-4 











1 








2 


-5 


-1 


1 


6 














1 





8 


6 


2 


-9 


1 

















1 



(5.4) 













Dual 












c 


W\ 


W-i 


w 3 


Wi 


w 5 


w$ 


fli 


V2 


v 3 


Va 





1 


3 


7 


5 


2 


8 














2 


4 


-1 


2 


3 


-5 


6 


-1 











13 


1 


3 


-4 


5 


-1 


2 





-1 








7 


3 


7 


1 


-2 


1 


-9 








-1 





6 


5 


-3 


8 


-4 


6 


1 











-1 



Carrying out a phase / step of the simplex algorithm on 5.4 gives us the dual tableau (see Table 5.5) 
in basic feasible form (the objective row has been omitted since it is not needed in the rest of the 
procedure). 

Now we construct an objective row for 5.3 from the basic feasible solution to the dual (5.5). This 
is shown in Table 5.6. 

After adjoining 5.6 to 5.3, four pivots of the simplex algorithm are required to minimize. The result 
is shown in Table 5.7. 

Since the objective function has not been reduced to zero (entry (1, 1) in 5.7), at least one dual 
cycle is required. We construct an objective row for the dual (5.5) from the solution to the primal (5.7). 
This is shown in Table 5.8. 

After adjoining 5.8 to 5.5, one pivot of the simplex algorithm is required to minimize. The result 
is given in Table 5.9. 

Since the value of the objective function in 5.9 is now zero, 5.9 must be the optimal dual tableau 
(except for the objective row, of course). Thus the previous primal tableau (5.7) must have been optimal. 
Checking the values of the actual objective functions using the solutions given in 5.7 and 5.9 for the 
primal and dual, respectively, we find that both are 13. 



MUTUAL PRIMAL-DUAL /LGORITHM 



205 



ACKNOWLEDGEMENTS 

The author wishes to thank Dr. Richard Young, Associate Professor of Economics, Rice University, 
and Mr. Arthur Ziffer of the U.S. Naval Research Laboratory for their many helpful suggestions toward 
the development of this algorithm. 

Table 5.5 



c 


w, 


W2 


w 3 


w 4 


W 5 


We 


V\ 


V2 


Va 


v 4 






1.2979 


1.0000 


0.0000 


0.8511 


0.0000 


0.0000 


0.3830 


-0.0164 


-0.0213 


-0.0426 


-0.0851 


2.2098 


0.0000 


0.0000 


-0.7328 


1.0000 


0.0000 


1.0040 


0.0247 


-0.2044 


0.0797 


-0.0267 


0.8179 


0.0000 


1.0000 


-0.4206 


0.0000 


0.0000 


-1.1544 


0.0351 


-0.0302 


-0.1069 


0.0421 


1.8006 


0.0000 


0.0000 


-0.0747 


0.0000 


1.0000 


-0.0604 


0.1227 


-0.1336 


0.0351 


-0.0925 



Table 5.6 



b 


Xl 


x 2 


X3 


Xi 


5l 


52 


5 3 


Si 


5 5 


5 6 


18.4018 


0.0000 


0.0000 


0.0000 


0.0000 


1.2979 


0.8179 


0.0000 


2.2098 


1.8006 


0.0000 



Table 5.7 



6 


Xi 


X 2 


x 3 


X4 


5l 


S2 


S3 


S4 


$5 


56 


5.4018 


0.0000 


0.0000 


-51.2537 


0.0000 


-3.4789 


-6.2958 


0.0000 


1.8733 


0.0000 


0.0000 


0.0000 


1.0000 


0.0000 


-3.4930 


0.0000 


-0.0423 


-0.4085 


0.0000 


0.2535 


0.0000 


0.0000 


.0.0000 


0.0000 


0.0000 


2.6338 


1.0000 


0.1972 


0.2394 


0.0000 


-0.1831 


0.0000 


0.0000 


11.0000 


0.0000 


0.0000 


2.1268 


0.0000 


-0.7606 


0.6479 


1.0000 


0.5634 


0.0000 


0.0000 


1.0000 


0.0000 


1.0000 


3.8028 


0.0000 


0.1831 


0.4366 


0.0000 


-0.0986 


0.0000 


0.0000 


3.0000 


0.0000 


0.0000 


-28.4648 


0.0000 


-1.2113 


-3.0423 


0.0000 


2.2676 


1.0000 


0.0000 


6.0000 


0.0000 


0.0000 


1.7183 


0.0000 


-0.3099 


1.3380 


0.0000 


-1.1408 


0.0000 


1.0000 



Table 5.8 



c 


Wi 


U>2 


W 3 


Wa 


M>5 


We 


V, 


V2 


va 


V4 


5.4018 


0.0000 


0.0000 


11.0000 


0.0000 


3.0000 


6.0000 


0.0000 


1.0000 


0.0000 


0.0000 











Table 5.9 












c 


W\ 


Ul2 


w 3 


W4 


(#5 


We 


Vl 


V2 


Va 


V4 


0.0000 


0.0000 


0.0000 


-11.0000 


0.0000 


-3.0000 


-6.0000 


0.0000 


-1.0000 


0.0000 


0.0000 


2.8593 


1.0000 


0.0000 


0.7863 


0.0000 


0.8672 


0.3306 


0.0000 


-0.1372 


-0.0122 


-0.1653 


1.8473 


0.0000 


0.0000 


-0.7178 


1.0000 


-0.2013 


1.0162 


0.0000 


-0.1775 


0.0726 


-0.0081 


0.3028 


0.0000 


1.0000 


-0.3992 


0.0000 


-0.2861 


-1.1371 


0.0000 


0.0080 


-0.1169 


0.0686 


14.6748 


0.0000 


0.0000 


-0.6088 


0.0000 


8.1500 


-0.4923 


1.0000 


-1.0888 


0.2861 


-0.7539 



REFERENCES 



[1] Balinski. M. L. and R. E. Gomory, "A Mutual Primal-Dual Simplex Method." in Recent Advances in Mathematical Program- 
ming (R. L. Graves and P. Wolfe, editors) (McGraw-Hill Book Co., New York, 1963), pp. 17-26. 

[2] Dantzig, G. B., "Composite Simplex-Dual Simplex Algorithm-I," The RAND Corp., RM-1274 (Apr; 1954). 

[3] Dantzig, G. B., L. Ford, and D. Fulkerson, "A Primal Dual Algorithm for Linear Programs," Paper 7 in "Linear Inequalities 
and Related Systems," Annals of Mathematical Studies (Harold Kuhn and Albert Tucker, editors) (Princeton University 
Press, Princeton, 1956), No. 38, pp. 171-181. 

[4] Dantzig, G. B., Linear Programming and Extensions (Princeton University Press, Princeton, 1963), Ch. 5 and 6. 



206 M - Y - HARRIS 

[5] Hadley, G., Linear Programming (Addison- Wesley, Reading, Mass., 1962), pp. 257-266. 

[6] Orchard-Hays, W., "A Composite-Simplex Algorithm-II," The RAND Corp., RM-1275 (May 1954). 

[7] Vajda, S., Mathematical Programming (Addison-Wesley, Reading, Mass., 1961), pp. 108-113. 



INTERCEPTION IN A NETWORK* 



R. D. Wollmer 

The RAND Corporation 
Santa Monica, California 



ABSTRACT 

This paper presents an algorithm for determining where to place intercepting units in 
order to maximize the probability of preventing an opposing force from proceeding from one 
particular node in an undirected network to another. 

The usual gaming assumptions are invoked; namely, the strategy for placing the units 
is known to the opponent and he will choose a path through the network which, based on 
this knowledge, maximizes his probability of successful traverse. As given quantities, the 
model requires a list of the arcs and nodes of the network, the number of intercepting units 
available to stop the opposing force, and the probabilities for stopping the opposition at the 
arcs and nodes as functions of the number of intercepting units placed there. From these 
quantities, the algorithm calculates the probabilities for placing the unit at the arcs and nodes 
when one intercepting unit is available, and the expected numbers of units to place at the 
arcs and nodes when multiple intercepting units are available. 

I. INTRODUCTIOIV 

The paper depicts a game between two opposing forces — an evader and an interceptor. Only the 
interceptor's problem will be dealt with. In this game, the evader attempts to proceed from one point 
to another in a network. The interceptor, who may possess one or more indivisible intercepting units, 
attempts to stop him by placing these units along arcs or nodes that he expects the evader to traverse. 
His problem is to place these so as to minimize the evader's probability of successful traverse. The 
evader's problem is, of course, to select a path that will maximize this probability. Both problems can 
be represented by a zero-sum two-person game. As the next section will show, however, the game matrix 
is very large and difficult to generate. By using a steepest ascent approach, a solution may be obtained 
much more easily for the interceptor. 

This paper presents an algorithm based on such an approach. As inputs, it requires a list of the 
arcs and nodes of the network, the number of intercepting units available to stop the opposing force, 
and the probabilities for stopping the opposition at the arcs and nodes as functions of the number 
(integer) of intercepting units placed there. From this information, the algorithm calculates the expected 
number of intercepting units to place at each arc and each node. For the case where only one inter- 
cepting unit is available, the expected values are probabilities of placement; for more than one intercept- 
ing unit, expectations are to be interpreted in the most obvious way. Specifically, if the expected number 
of units to place at a particular location is 4-1/4, one would always place four there and would place 
a fifth one there one-fourth of the time. 

The situation depicted by this algorithm is particularly relevant in applying interdiction to infiltra- 
tion and counterinsurgency, where the counterinsurgent or infiltrating force would play the role of the 



*This research is supported by the United States Air Force under Project RAND. Views or conclusions contained in this 
study should not be interpreted as representing the official opinion or policy of the United States Air Force. 

207 



208 R- D. WOLLMER 

evader, while the interdiction forces would play that of the interceptor. It also applies to police networks 
where a wanted fugitive (evader) is attempting to escape from a city while police (intercepting) units 
are attempting to close off his avenues of escape. 

The algorithm to be presented is a steepest ascent approach and yields solutions that are always 
optimal in the gaming sense (i.e., the evader's best chance of successful traverse is as small as possible) 
for the interceptor when he has but one intercepting unit at his disposal. For the case of multiple inter- 
cepting units, the solutions are not always optimal; thus for the general case the algorithm is an 
approximation scheme. 

II. INTERCEPTOR'S PROBLEM FOR ONE INTERCEPTING UNIT 

The network is characterized by sets of elements called arcs and nodes. Nodes are points or 
junctions, and arcs are undirected line segments joining nodes. An arc that joins node i to node j may 
be disignated by the symbol (i,j) or (J, i). The node from which the evader starts is called the source; 
the one he attempts to reach is called the sink. 

It is assumed that the evader and interceptor are each composed of one indivisible unit. Later on 
this restriction will be relaxed for the interceptor. The evader is attempting to choose a path from 
source to sink in such a way that his probability of successful traverse is maximized. The interceptor 
wishes to place his intercepting unit so that this probability is minimized. As is customary in game 
theory, it is assumed that each side's strategy will become known to the other and consequently each 
will need to adopt a mixed strategy. The inadequacy of pure strategies for the interceptor may be dem- 
onstrated by any network whose source and sink are relatively invulnerable to attack and any particular 
intermediate arc or node may be bypassed. Placing the intercepting unit at the source or sink would be 
relatively ineffective while placing it anywhere else would be completely ineffective as the evader 
would choose a path that bypassed the unit. 

Quantities defining the local effectiveness of placing the intercepting unit at arcs and nodes, 
and quantities expressing the interceptor's strategy are as follows: 

p(i) = probability that the evader will be stopped at node i, given that node i is on his chosen path 
and the interceptor places his intercepting unit there. 

p(i, j) = probability that the evader will be stopped at arc (i, j) given that arc (i, j) is on his 
chosen path and the interceptor places his intercepting unit there. 

77 (i) — probability that the interceptor places his intercepting unit at node i. 

7r(i,y) = probability that the interceptor places his intercepting unit at arc (i,j). 

The probability that the evader can successfully cross a particular node, t, is 1 — n(i)p(i) and 
that of crossing a particular arc (i,j) is 1 — Tr(i,j)p{iJ). Assuming these probabilities are independent* 
and the evader attempts to reach the sink by the path i,, . . . i n , his probability of successful traverse is 

n n-1 

(1) ^ = E[ [ l-1T (ir)p(ir)] ][ [1— ir(ir, ir+l)p(ir, U+l)]- 

r=l r=\ 

This quantity, K, will be referred to as the value of the path ii, . . ., l». 

Given the p{i), p(i,j), tt(i), and v(i,j), the evader's problem is to find a source-sink path of 
maximum value. The interceptor's problem, given the p(i) and p(i, j), is to choose 7r(t) 2*0 and 



*This assumption is in fact not true. To illustrate this, suppose ir(i) = l/2 and n(j) = 1/2. Independence would require 
that the probability of the intercepting unit being at nodes i and j at the same time be 1/4 while in reality it is zero. However, 
the next section will show that this assumption leads to no inaccuracies for the one-intercepting-unit case. 



INTERCEPTION IN A NETWORK ._ 209 

7r(i, j) 3 s such that the maximum value of all source-sink paths is minimized, subject to the con- 
straint %Tr(i) +2,ir(i,j) *£ 1. Note that if all p(i) and p(i,j) are strictly positive, any optimal strategy 
will be such that strict equality holds in this relationship. To see this, suppose the intercepting unit 
is not placed at any arc or node with probability p (i.e., 2/7r(t) + 2/7r(i, j) = 1 — p). Then increasing 
all ir(i) and Tr(i,j) by equal amounts totaling p yields an improved solution. 

III. A STEEPEST ASCENT APPROACH FOR ONE INTERCEPTING UNIT 

The problem of finding the n{i) and ir(i,j) can in theory be solved by linear programming. The 
linear program solution would include an optimal strategy for the evader as well as the interceptor. 
The value of the game would be equal to both the evader's maximum guaranteed probability of 
successful traverse, and one minus the interceptor's guaranteed probability of stopping the infiltrator. 
However, the game matrix for the program would require a column for every possible source-sink 
path and a row for each of the interceptor's pure strategies. The number of source-sink paths is exceed- 
ingly large even for relatively small networks and in later sections, when multiple intercepting units 
are allowed, the same may be said of the number of pure interceptor strategies. 

The algorithm of this section is a steepest ascent approach which avoids these enumeration 
problems; however, while it yields an optimal strategy for the interceptor and a value for the game, it 
does not yield an optimal strategy for the evader. 

Initially, all 7r(i) and ir{i,j) are assigned values of zero. Then they are increased by small amounts 
A7r(i) 2=0 and A-7r(i, j) 2=0. The proportions for the Air(i) and ^rr(i,j) are such that the additional 
intercepting unit placement probability, 2A7r(£) +XA7r(i, _/), divided by the decrease in maximum 
path value, is minimized as these two quantities tend toward zero. This is equivalent to maximizing 
the decrease in maximum path value per unit of additional intercepting unit probability allocation. 
The 7r(t) and ir{i, j) are then increased again in the same manner until they sum up to one. 

For a specific path, t'i, . . ., i n , increasing Tr(i r ) by A7r(i r ) replaces the factor [1 — ir(i r )p(i r )] 
in Eq. (1) by [1 — (rr(i r ) + A7r(£ r ))p(£ r )J and reduces the value of the path from K to K — &K, where 

1 l-ir(»r)p(»r) J" 

Solving for the quotient of the additional intercepting unit probability allocation and the decrease in 
path value, the following expression is obtained: 



(6a) AK K 



ir(ir) 



P(lr) 



If intercepting unit placement probability were increased at an arc instead of a node, the same 
argument would yield an expression for the additional intercepting unit probability placement divided 
by the decrease in path value which is identical to the preceding one, with arc probabilities substituted 
for node probabilities. Specifically, for arc (i r , ir+i), this would be 



(3b) **&' ir+l) ~ 



7T(i r , Ir+l) 



A£ Klp(i r ,i r+1 ) 

Note that in order to obtain any decrease in maximum path value, the A7r(£) and kir{i,j) must be 
strictly greater than zero along a subset, C, of arcs and nodes that interesects all paths of maximum 



210 R- D. WOLLMER 

value. An efficient allocation would require that all members of C remain on maximum value paths, 
for otherwise, lessening the probability increase at an arc or node not on a maximum value path by a 
small amount would decrease SA7r(i) -+- 2A7r(i, j) , while not affecting the maximum path value. Then, 
of course, this small amount of probability could be redistributed among arcs and nodes on maximum 
value paths to get a strict improvement. Thus, the Air(j) and k-rr(i,j) assigned to the chosen C must be 
such that the A£'s are equal for each of its members. In other words, they must be proportional to 

77(0 



(4) p(i) 

b) , -7T(i,j). 

p(l,J) 

The constant — , where K is the maximum path value, was dropped in Expression (4), of course. The 

total increase in intercepting unit placement probability per unit decrease in the value of a maximum 
probability path would be: 






-n(i,j) 



p(ij) 



Thus, the problem is reduced to finding a subset of arcs and nodes intersecting all maximum prob- 
ability paths that minimizes Expression (5). 

Initially, all paths have value one since all ir(i) and 7r{i, j) are zero. Thus C is required to be a 
set blocking all source-sink paths. If the nodes and arcs are assigned capacities equal to the quantities 
of Expression (4a) and (4b), the problem of finding C is reduced to one of finding a minimum cut.* 
However, this is equal to the value of the maximum source-sink flow and can be found by the maximum 
flow algorithm of Ref. [3]. 

Note that once additional intercepting unit probability is allocated to the minimum cut, C, the 
Capacity of each arc and node in C decreases by its additional probability increase while the capacity 
of any other arc or node remains the same. Thus, the value of C will be reduced by SA7r(i) +SA7r(i, j) 
while all other cuts will decrease by amounts that do not exceed this. Hence, C remains minimum and 
the entire unit of probability may be allocated to C. 

The 7r(i) and ir{i, j) may therefore be solved for as follows: 

1. Assign all nodes capacities of llp(i) and all arcs capacities of l/p(i, j) (since all 7r(i) and 
7r(i, j) are zero). 

2. Maximize flow from source to sink. 

3. Let C be the minimum cut set and ^ its value. Set ir(i) =p(i)/V for ieC and Tr(i,j) =p(i,j)IV 
for (i,j)eC. Set all other 7r(i) and Tr(i,j) equal to zero. 

Note that since only one cut is obtained, each path of maximum value has positive intercepting 
unit probability allocated to only one of its arcs or nodes, and both its value and traverse probability 
are equal to the probability associated with that particular arc or node. Hence, the independence 
assumption does not lead to inaccuracies for the one-intercepting-unit case. Finally, the maximum path 

value is 1 — -p, where V is the value of C found in step 3. 



*Essentially, a cut set is a set of arcs and nodes blocking all source-sink paths. Its value is the sum of the capacities of its 
arcs and nodes, and the minimum value of all cuts is equal to the maximum value of the source-sink flow. See Ref. [3 J. 



INTERCEPTION IN A NETWORK 211 

IV. MULTIPLE INTERCEPTING UNITS 

The steepest ascent algorithm of the last section may be modified to handle the case where there 
is more than one intercepting unit; however, the interceptor strategy obtained in some cases may not 
be optimal and the game value obtained may be approximate rather than exact. This will be discussed 
more fully later on. The quantities measuring the local effectiveness of placing intercepting units at 
arcs and nodes and those defining the interceptor's strategy must be redefined as follows: 

p(i)k — probability that the evader will be stopped at node i, given that node i is on his chosen 
path and the interceptor places k intercepting units there. 

pi.ii J )k — probability that the evader will be stopped at arc (i,j), given that arc (i,j) is on his chosen 
path and the interceptor places k intercepting units there. 

tt{i) = expected number of intercepting units the interceptor places at node i. 

ir(i,j) — expected number of intercepting units the interceptor places at arc (i,j), 
where p(i)o andp(i,j)o are both identically zero. Also let 

7r(i')=7r(;) — [«■(*)] 
ir{i,j) =ir(i, ;') — [n(i, j) ] . 

Note that 7r(i) and 7r(i,j) are merely the fractional parts of 7r(i) and n(i,j). It will be assumed 
p(i)k and p(i,j)k are increasing in k and 

(6) a) p(i)k+i—p(i)k^p{i)k—p(i)k-i 

b) p(i, j)k+i—p(i, j)k < p(i, j)k~ p{i, j)k-u k > 1. 

The assumption of Expression (6) essentially states that if an additional intercepting unit is placed 
at any arc or node, it cannot be more effective in stopping the evader than any previous unit placed 
there. It is often satisfied in practice as the effects of overkill become more pronounced as the number 
of units at a location increase. The simplest example of this is the case where each unit placed at node i 
has independent probability p of stopping the evader. For this situation p(i)fc— 1 ~~ (1 — p) k and 
p(i)k+i~ p(i)k = p(l— p) k which is certainly decreasing in k, thus satisfy (6). 

The values of 7r(i) and n(i,j) of course are not sufficient in themselves to define the interceptor's 
strategy as any non-zero value of 7r(i) or ir(i,j) may often be realized in several ways. However, the 
assumption of Expression (6) assures that for a given value of -n ( i ) , the probability of stopping the evader 
at node i should he attempt to cross it is minimized by allocating [7r(t)] + 1 units at i with probability 
7r(i), and [7r(i)] units with probability 1 — Tr(i). In other words, if 7r(i) = 3-1/2, tt{i) would best be 
realized by assigning three forces at i half the time and four half the time, as opposed to such a policy 
as allocating two forces half the time and five half the time. This is shown in the theorem to follow. First, 
define P{k, i} as the probability that k units are placed at node i. Note that 2,kP{k, i} = ir(i) and the 
probability of thwarting an evader's attempt to cross node i is 1 — 2,P{k, i}p(i)k. Minimizing this prob- 
ability is equivalent to maximizing 2,P{k, i}p(i)k. 

THEOREM: The quantity (t')2P{/:, i}p(i) k is maximized, subject to the constraint (ii)^,kP{k, i} 
= ir(i) when P{[ir(i)] + 1, i}=--n-(0 ,<P{[n(i)], i} = l-n(i), and all other P{k, i} = 0. 

PROOF: Let h = max {k/P{k, i}>0}, /= min {k/P{k, i} > 0}, and d=h-l. Let P{k, i} = P{k, i} 
be a set of values of smallest d which maximizes (i) subject to (ii). Suppose the theorem is false. 
Then h-l = d^2. Let m = min{P{h, i}, P{1, i}}. Set P{h, i} = P{h, i} - m and P{h-1, i} 
= P{h — l, i} + m. Then set P{1, i} = P{l, i} — m and increase P{/-f 1, i} by m. (If h — l = 2, then 
P{l+l} = P{l+l} + 2m.) Either P{h, i} or P{1, i} is now zero, decreasing d by at least one unit. 



212 R- D. WOLLMER 

The increase in (/) is equal to m(— p(i)i + p(i)i+\ — p(i)ft + p(0/i-i)- Since (ii) is still satisfied, this 
is strictly negative from our choice of P{k, i}. However, since h > / + 1, 

P(i)h — p(i)h-i =S p(i)i+i—p(i)i 
- p(i)i + p(i)i+i— p(i)h + p(i)h-i 2= 

for a contradiction. Q.E.D. 

The same results, of course, hold for the Tr(i, j) of the arcs. The significance of this is that the 
search for an optimal interceptor strategy may be confined to strategies of the above type. When this 
is done, the 77(1) and Tr(i, j) are sufficient to specify the interceptor's strategy. 

As in Eq. (1) of Sec. II, if independence is assumed, the value of a path is still the product of the 
probabilities associated with its nodes and arcs; however, the probability associated with node i is now 

7T(j) (1— p(i)[ir(l)]+l) + (1-17(0) (1— P(i)[n(l)]) 

or 

^~p(i)[Mi)] — rr(i)(p(i)i7T(i)]+i—p(i)[Tr{i)]), 
and that with arc (i, j) is 

l^P(^y)[7r(i,j)] — 7r(l, j)(p(i, j)[7r(i,j)] + l— P(i, j)[n(i,j)]). 

If ti, . . ., in is a path of value K, and 7r(i r ) is increased by A7r(i r ), reducing the path value to 
K — &K, the expressions for K and K — &K, as before, differ only in the factor for node i r . Thus, for 
sufficiently small A7r(i r ), the expression for the quotient of the intercepting unit increase and path 
decrease is 

A7T(l r ) 1 f 1—P{ir)[n(i r )) 1 

-TT{l r ) \ , 



AX K [p(ir)[n(i r )]+l— P(ir)[7r(i r )] 

and decreasing path value by increasing the number of intercepting units at an arc instead of a node 
yields the analogous expression, 

A7T(ir, ir+l) If 1 — P(ir, ir+l)[n(i r ,i r+l )) 



A£ 



If 1 — PKlr, lr + l)[ir(i r ,i r+l )] ... 1 

= -p\ — — : : T. — : ^ ?rUr, ir+\) \ 

A I p{l r , lr+l)ln(i r ,i r + 1 )] + l—PKlr, lr+l) [n(i r , i r + 1 )] J 



As in Sec. Ill, intercepting units must be allocated along all arcs of some cut set, and the problem 
is reduced to finding a minimum cut set, which in turn reduces to one of finding a maximum source- 
sink flow. The node and arc capacities are now 

(7) a) C(l)=—7-^. 7TT 77 (I) 

P(l)[77(t)]+1-P(l)[7r(i)] 

b) C(l,J)= ,. ., /• -x TT\l,J)- 

Note that the c(i) and c{i, j) are no longer decreasing functions of 77(1) and 77(1,7) except over 
intervals whose endpoints possess the same integer part. Specifically, c(i) is a decreasing function of 



INTERCEPTION IN A NETWORK 213 

7r(t) for < -n-(i) < 1, 1 s£ 7r(i) < 2, etc. However, c(i) is discontinuous and actually increases at all 
integer points for which (6a) holds with strict inequality. This follows since 

lim l—p(t) [„(!)] . ,., _ 1— ■p(Q[ir(i)]+i < 1— p(t)[7r(»)]+i 

ff(i) ^ 1 p(i)[n(i)]+l— P(i)[ir(i)] P(i)[n(i)]+1—P(i)[n(i)] P (0 [n(i)]+2 ~ P (0 [ir(i)]+l 



Similar results hold for c{i, j). Thus, when a minimum cut set is found and additional intercepting 
units are allocated along its arcs and nodes, that cut remains minimum provided the expected number 
of intercepting units on these arcs and nodes do not reach or exceed their next highest integer values. 
This requires that the steepest ascent approach be modified as follows. 

Starting with all 7r(t) and ir(i, j) equal to zero, capacities are assigned to the arcs and nodes of 
the network which are equal to the quantities defined by Expression (7). Flow is maximized from 
source to sink and the 7r(t) and Tr(i,j) corresponding to the arcs and nodes of the minimum cut set are 
increased by amounts A7r(i) and A7r(i, j) which are proportional to their capacities. The constant of 
proportionality, M, is the smallest possible such constant that will either increase some 7r(t) or Tr(i,j) 
to the next highest integer or will increase %Tr(i) + ?.TT{i,j) to n. If the latter does not happen, the new 
values of ir(i) and rr(i,j) are used to calculate new capacities and flow is maximized again to obtain 
a new cut. The process is repeated until %7r(i) +£77- (1,7) = n. 

Specifically, the algorithm for the general case of n intercepting units is as follows: 

1. Set all 7r(i) = 0and all ir{i,j)=0. 

2. Assign the nodes and arcs of the network capacities of c(i) and c(i,j) respectively, where c(i) 
and c(i, j) are as defined by Expression (7). 

3. Maximize flow from source to sink and let C be the nodes and arcs of the minimum cut set found. 

A n Tit {1 _ tHO} 

4. Compute M x = mm - 1 ... 

i€C c(i) 

1VI2 — mm — — 

(i,j*C C(l,j) 

5. Let M=min (M u M 2 , M 3 ). Set bn(i)=Mc(i) if ieC and &ir(i,j)=Mc(iJ) if (i,j)eC. Set all 
other A7r(t) and Air(i,j) equal to zero. 

6. Set ir(i) = ir{i) + Att({) and ir(ij) = n(i,j) + bir(i,j) for all i and (i,j). If %ir(i) = £w<i,» < n, 
go back to step 2. Otherwise terminate. 

Step 1 is simply an initialization step. Steps 2 and 3 determine which nodes and arcs should have 
their expected number of intercepting units increased and the proportions of the increase. Specifically 
the expected number of intercepting units will be increased along the nodes and arcs of the minimum 
cut set found in proportion to their capacities. Steps 4 through 6 determine the proportionality constant 
M, makes the appropriate increases, and makes a termination test. Note that if M = M\, some -rr{i) 
will be increased by 1 — Tr(i), and since 77(1) is the fractional part of 7r(i), this will increase 7r(i) to an 
integer value. Similarly if M=?M 2 , some v(i, j) will be increased to an integer value. If M = M 3 , total 
expected intercepting units will be increased by n-XTr(i) —%TT(i,j) to n. In the latter case all n inter- 
cepting units have been allocated and the procedure terminates. 



214 



R. D. WOLLMER 



At termination the ir{i) and rr{i, j) represent the expected number of intercepting units to place 
at node i and arc (i, j), respectively. As mentioned earlier, the strategy will be to place [77(1)] + 1 
forces at node i with probability ir(i) and [77(1)] forces with probability 1 — Tr(i). 

The maximum path value can be kept track of during the course of the algorithm. Specifically, 
if K is the maximum path value at the beginning of an iteration, then the maximum path value at the 
end of the iteration is K(l — M), where M is that found in Step 5. Of course, K= 1 at the beginning of 
the first iteration. This is easily proven by letting V be the value of C, the minimum cut set. Then 
substitute the generalized capacities (Expression 7a— b) for the specialized ones (Expression 4a— b) 
in Expressions (3) and (5). These replacements for Expressions (3) and (5) will have the same meaning 
in the multiple intercepting unit situation as Expressions (3) and (5) did for the single intercepting 
unit situation. The replacement for (5), which represents the total increase in intercepting unit place- 
ment probability per unit decrease in the value of a maximum probability path, is: 

^Air(i) +SA-n-(t, /) = 1 y 
AK K 

However, 2A7r(i) +£A7r(i, j) =MV where M is as defined in Step 5 of the multiple intercepting unit 
algorithm. Thus 

MV = \y 

AK K 



or 



K-AK = K(l-M). 



Note that when there is more than one intercepting unit to be assigned, the strategy obtained is 
not necessarily optimal. To illustrate this, consider the network of Fig. 1. 




Figure 1. A simple network with node 1 the source and node 6 the sink. 

Suppose p(l, 2)i = 5/8, p(l, 2)2=1, p(2) fc = l — (1/4)*, and all other p(i)k and p(i, j) k equal 
1 — (1/2)*. If there were two intercepting units available, placing both at arc (1, 2) with probability 
one will make the evaders successful probability of traverse equal to zero; however, the algorithm, as 
may be easily verified, will place both units at node 2 with probability one giving the evader a successful 
traverse probability of 1/16. For one unit, however, the algorithm will optimally place it at node 2 
giving the evader a successful traverse probability of 1/4. 



N 




4 8 

FIGURE 2. Sample problem network with node 1 as source and node 9 as sink. 



INTERCEPTION IN A NETWORK 215 

V. EXAMPLE 

For illustrative purposes, the algorithm is used to find the optimal placement of forces for the 
network of Fig. 2, with the results of each iteration. Ten intercepting units are available, and all arcs 
and nodes have identical probabilities associated with them. Thus, for each node, i, and each arc 
(i, /), we have: 



p( 
p( 
p( 
p( 
p( 
p( 
P( 
P( 
P( 



) , = 0.6666 p(i,j)i = 0.6666 

) 2 = 0.8332 p(i,j) 2 = 0.8332 

) 3 = 0.8748 p(i,» 3 = 0.8748 

) 4 = 0.8852 p(i,y) 4 = 0.8852 

) 5 = 0.8878 p(i,y) 5 = 0.8878 

)« = 0.8884 p(i,y) fi = 0.8884 

) 7 = 0.8886 p(i, j) 7 =0.8886 

) H = 0.8887 p(i, y)« = 0.8887 

) 8 = 0.8888 p(i,» » = 0.8888 



p(0 ,0 = 0.8889 p(i,j) ,0 = 0.8889 

The problem terminates in eight iterations. The results of these iterations follow: 

Iteration 1: All arcs and nodes have capacities of 1.5. A maximum flow is 1 unit along path 1,2,5,6,9. 
A minimum cut set is node 1 and its value is 1.5. Increasing intercepting units at node 1 to 1 gives a 
maximum path value of 0.33333. 

Iteration 2: The capacity of node 1 is increased to 2.0; all other capacities are unchanged from 
iteration 1. A maximal flow is 1.5 units along path 1,2,5,6,9. A minimum cut set is node 5 and its value 
is 1.5. Increasing intercepting units at node 5 to 1 gives a maximum path value of 0.11111. 

Iteration 3: The capacity of node 5 is increased to 2.0; all others are unchanged. A maximal flow 
is 1.5 units along path 1,2,5,6,9. A minimum cut set is node 9 and its value is 1.5. Increasing intercepting 
units at node 9 to 1 gives a maximum path value of 0.03704. 

Iteration 4: The capacity of node 9 is increased to 2.0; all others are unchanged. A maximal flow 
is 1.5 units along path 1,2,5,6,9 and 0.5 unit along path 1,3,5,8,7,9. A minimum cut set is node 1 and 
its value is 2.0. Increasing intercepting units at node 1 to 2 gives a maximum path value of 0.01853. 

Iteration 5: The capacity of node 1 is increased to 4.0; all others are unchanged. A maximum flow 
is 1.5 units along path 1,2,5,6,9 and 0.5 unit along path 1,3,5,8,7,9. A minimum cut set is node 5 and 
its value is 2.0. Increasing intercepting units at node 5 to 2 gives a maximum path value of 0.00927. 

Iteration 6: The capacity of node 5 is increased to 4.0; all others are unchanged. A maximal flow 
is 1.5 units along path 1,2,5,6,9 and 0.5 unit on path 1,3,5,8,7,9. A minimum cut set is node 9 and its 
value is 2.0. Increasing intercepting units at node 9 to 2 gives a maximum path value of 0.00464. 

Iteration 7: The capacity of node 9 increased to 4.0; all others are unchanged. A maximal flow is 
1.5 units along path 1,2,5,6,9 and 1.5 units along path 1,3,5,8,7,9. A minimum cut set is arcs (5,6) and 
(5,8) and its value is 3.0. Increasing intercepting units to 1 at arcs (5,6) and (5,8) gives a maximum path 
value of 0.00155. 

Iteration 8: The capacities of arcs (5,6) and (5,8) are increased to 2.0; all others are unchanged. 
A maximal flow is 1.5 units along path 1,2,5,6,9 and 1.5 units along path 1,3,5,8,7,9. A minimum cut set 
is nodes 6 and 8 and its value is 3.0. Increasing intercepting units to 1 at nodes 6 and 8 gives a maximum 
path value of 0.00052. Algorithm terminates as all 10 units have been assigned. 



216 R D - WOLLMER 

The final assignment is 2.0 intercepting units each at nodes 1,5, and 9; and 1.0 intercepting unit 
each at nodes 6 and 8 and arcs (5,6) and (5,8). 

REFERENCES 

[1] Durbin, E. P., "An Interdiction Model of Highway Transportation," The RAND Corporation, RM-4945-PR (May 1966). 

[2] Ford, L. R., Jr. and D. R. Fulkerson, "Flows In Networks," The RAND Corporation, R-375 (Dec. 1960). 

[3] Wollmer, R. D., "Maximizing Flow Through A Network With Node and Arc Capacities," Transportation Science 2, 213-232 

(1968). 
[4] Wollmer, R. D., "Removing Arcs From A Network," Operations Research 12, 934-940 (1964). 
[5] Wollmer, R. D., "Stochastic Sensitivity Analysis of Maximum Flow and Shortest Route Networks," Management Science 

9, 551-564 (1968). 
[6] Wollmer, R. D., "Algorithms For Targeting Strikes in a Lines of Communication (LOC) Network", Operations Research, 18, 

497-515 (1970). 



THE FIXED CHARGE PROBLEM* 



David I. Steinberg 



Washington University 
St. Louis, Missouri 



ABSTRACT 

The fixed charge problem is a nonlinear programming problem of practical interest in 
business and industry. Yet, until now no computationally feasible exact method of solution 
for large problems had been developed. 

In this paper an exact algorithm is presented which is computationally feasible for large 
problems. The algorithm is based upon a branch and bound approach, with the additional 
feature that the amount of computer storage required remains constant throughout (for a 
problem of any given size). 

Also presented are three suboptimal heuristic algorithms which are of interest because, 
although they do not guarantee that the true optimal solution will be found, they usually 
yield very good solutions and are extremely rapid techniques. 

Computational results are described for several of the heuristic methods and for the branch 
and bound algorithm. 

1. INTRODUCTION 

The fixed charge problem is perhaps one of the most intriguing problems of mathematical pro- 
gramming. Its structure is almost identical with that of a linear programming problem; yet, the one 
difference — the existence of the fixed charges in its objective function — has prevented the development 
of any extensive theory concerning its solution. 

There are two known exact methods of solution: total enumeration of the extreme points of the 
convex set of feasible solutions, and a mixed integer-continuous variable formulation. However, total 
enumeration is impractical for large problems, and, at this time, methods for solving mixed integer- 
continuous variable problems involving a large number of integer variables are not effective enough 
to be of much use, either. Therefore, heuristic methods for obtaining a "good" (although not necessarily 
optimal) solution have been developed. Some of these methods will be discussed in Section 2. 

Also presented is an exact algorithm, based on a branch and bound approach, which has been 
developed by the author, along with several new heuristic methods. Computational results are given 
in Section 4. 

1.1 Formulation of the Fixed Charge Problem 

The fixed charge problem is a nonlinear programming problem of the form: 

(1) Minimize z = cx + dy 

subject to: Ax — b, 

(2) x ^ 0, 



This research was performed under AEC Research Contract No. AT(11-1)-1493. Submitted originally as Report No. 
COO-1493-15. 

217 



218 D. I. STEINBERG 

where A is an m X n matrix of real numbers, 

b is an m-vector in an m-dimensional space, E m , 

c and d are n-vectors in E", 

x is an rc-vector of unknowns, and 

y is an n-vector such that 

fO, if*, = 

^ ^={i, if,;> - 

It is assumed that all dj ^ 0. The quantities dj are the fixed charges which are incurred if the corre- 
sponding Xj are (strictly) positive. If all dj = 0, then (1) and (2) reduce to a linear programming problem. 
It is known [7] that the optimal solution to the fixed charge problem occurs at an extreme point 
of the convex set of feasible solutions (2). Since the total number of extreme points is finite, the problem 
could be solved by evaluating the objective function (1) at every extreme point and taking the minimum 
of these quantities. However, as is the case with linear programming, the number of extreme points 
increases combinatorially with m and n. Hence, it is desirable to develop a method that requires that 
only a subset of the extreme points be examined in order to determine the optimal solution. Two such 
methods are the mixed integer-continuous variable integer programming formulation and the author's 
branch and bound algorithm. 

1.2 A Mixed Integer-Continuous Variable Approach 

Hadley [7] presents the following formulation of the fixed charge problem and suggests that it be 
solved by use of Gomory's mixed integer-continuous variable algorithm. 

Minimize z = cx + dy, 

subject to: Ax = b 

Xj-Ujyj^O, 7=1, . . ., n 
0<ft<l, 7-1 n 

yj integers 
x^O, 

where the Uj represent known upper bounds for the xj. 

If A ismXn, then the above formulation has m + 2n constraints. It should be noted, however, that 
the upper bound constraints on the yj need not be expressed explicitly, since no optimal solution would 
have yj > 1. Thus, the actual basis size required to use Gomory's algorithm is m + n. However, compu- 
tational experience with Gomory's algorithm has been very disappointing, even for some fairly small 
problems [7]; therefore, this approach is not widely used to solve fixed charge problems. 

2. SOME HEURISTIC METHODS 

The fact that the optimal solution to the fixed charge problem lies at an extreme point of the 
convex set (2) has led to the development of some heuristic approaches which use adjacent extreme 
point techniques to obtain good solutions. All of the methods discussed below are basically adjacent 
extreme point methods. 

There are also several heuristic methods for solving the fixed charge transportation problem. 
These are presented in [2] and [3], but shall not be discussed here. 



FIXED CHARGE PROBLEM . 219 

The following notations are offered prior to presenting the heuristic methods: 

Xb~ vector of basis variables, 

c B = vector of prices for variables in the basis, 

aj=yth column of A, 

B= basis matrix, 

yj — B~ l a,j, and 

From the theory of linear programming [6, 7] it is known that once a vector a,- is chosen to enter 
the basis, the vector a B r which leaves the basis is determined by 

XBr 

(4) — = Min { XBklykj , for y kj > 0} . 

Yrj k 

If all dk = (i.e., the problem is all linear) the greatest decrease in the objective function at each 
iteration is obtained if the vector to enter the basis, a,, is chosen so that 



X t 



(cj — Zj) =Min | — {c k — z k ), for y rk > [ 
k [yrk J 



When this expression is nonnegative the optimal solution has been reached. For the fixed charge 
problem (one or more dj ^ 0) criterion (5) must be modified to account for the fixed charges; the greatest 
decrease in z is obtained if a, is chosen so that 

(6) ej = Min {e fc , for y rk > 0}, 

k 

where 

(7) e k = — (ck~ z k ) + (d Br — d k ) 

Vrk 

Equation (7) must be modified slightly in the presence of degeneracy, in an obvious way. (It should 
be noted that the variable which leaves the basis, x Br , will, in general, be different for each,/). Because 
the objective function is concave, it is not true that when all e* 2s a global minimum has been reached. 
This stage may merely be a local minimum. That is, no adjacent extreme point yields an improved 
value of r, but it is still possible that some other extreme point of the convex set will yield a better value 
of 2. The methods discussed below differ in their approach to the problem of trying to obtain a better 
solution once a local minimum has been reached. The first three methods are those of the author. 

2.1 Heuristic One 

1. By using a Simplex Algorithm, vectors are removed from and entered into the basis according 
to (4) and (6), respectively, until all e,- 2= 0. Let the current solution be denoted by xo, with corresponding 
value of the objective function zo. Store xq. Define two parameters: 

ao = maximum number of consecutive iterations allowed to find a new extreme point with a value 
of 2 smaller than zo. 

/3o = maximum number of consecutive iterations allowed in which the objective function does not 
decrease. 

Also, define two "counters-," ai and /3i, which indicate the current number of consecutive itera- 
tions for which zo has not been improved, and in which the objective function does not decrease, respec- 
tively. Set ai = /3i = 0. 



220 D. I. STEINBERG 

2. Find the e, which yields the smallest increase in z and insert a, into the basis. Set ai = /3i = l. 
Proceed to Step 3. 

3. Find the e, which yields the largest decrease in z and insert a, into the basis. Proceed to Step 4. 

4. Compare the new z, z\, with zo: 

a) If Zi < Zo, return to Step 1. 

b) If z\ = Zo, but Xi 9^ Xo, set ai = ai + 1. If ai =£ ao return to Step 3. Otherwise terminate. 

c) If Z\ = Zo and *i = x , set «i = ai + 1 and /3i = /3i + 1. If ai =£ « and /3i =£ /3 go to Step 5. Otherwise 
terminate. 

d) If z\ > zo, set ai = ai + 1. If ai =£ ao return to Step 3. Otherwise terminate. 

5. Perform /3i consecutive interations in each of which the vector which yields the largest increase 
in z is inserted into the basis. Return to Step 3. 

This algorithm attempts to find better solutions by moving away from the current local minimum 
in as few iterations as possible. The first such attempt (Step 2) actually increases the objective function 
as little as possible; if this attempt fails, the objective function is increased as much as possible (Step 
5) with the hope that upon once again decreasing (Step 3) a better local minimum will be found. The 
efficiency of this method, as well as its accuracy, depends on the parameters ao and /3o. In general, ao 
should probably be somewhat larger than /3o. Since computational experience [4] has shown that the 
average number of iterations required to solve a linear programming problem with m constraints is 
about 3/2m, a reasonable range for a might be a e[m/2, m]. As the computational results of Section 4 
indicate, much smaller values may also be used fairly successfully. Of course, economic considerations 
may also play a part in the determination of ao and j8o. 

2.2 Heuristic Two 

1. Same as in previous heuristic, except that the entire Simplex tableau corresponding to xo is 
saved. Define a as before, set /3 = n — m, ai = l, /3i=l, and proceed to Step 2 (assuming again that 
all ej 3= at xo). 

2. Beginning with the tableau corresponding to xo, find the e, which yields the 03i)th smallest 
increase in z and insert a, into the basis. Set /3i = /3i + 1, and if /3i =S (J3o + 1) proceed to Step 3. Other- 
wise terminate. 

3. Find the ej which yields the largest decrease in z and insert a, into the basis. Set ai = ai + l 
and if ai =£ ao, proceed to Step 4. If ai > a , proceed to Step 2. 

4. Compare z\ with zo: 

a) If z\ > zo, return to Step 3. 

b) If z\ — zo, but jci t^ xo, return to Step 3. 

c) If Zi = zo and x\~ xo, set ai = 1 and return to Step 2. 

d) If z\ < zo, return to Step 1. 

2.3 Heuristic Three 

Heuristic Three employs a Monte Carlo approach to find an extreme point with a better value of 
z than xo (where again at xo all ej 3= 0). In this method, extreme points within a certain "distance" 
from Xo are randomly examined. "Distance" here shall mean the number of basis changes required to 
move from one extreme point to another. Let 8 = maximum distance from xo allowed and let ao again 
be the maximum number of iterations allowed to find an improvement in z. 

The procedure is as follows: 

1. Insert a.j into the basis, where ej is the maximum decrease in z. Continue in this manner until 



FIXED CHARGE PROBLEM 221 

all e, S 5 0. Denote the current extreme point by x . Set a.\ =0 and proceed to Step 2. 

2. Randomly select (according to the procedure described below) a vector within a distance 
8 of Jto. Denote this vector by X\. Set a t = a t + 1 and proceed to Step 3. 

3. Compare z\ with zo: 

a) If z\ > zo and «i 3= a , return to Step 2. If z x > zo and «i > a , terminate. 

b) If z\ ^ zo, return to Step 1, replacing xo and zo by *i and zi, respectively. 

In order to complete the description of the above algorithm, the procedure for the random selection 
of extreme points must be described. For simplicity, the vectors selected to enter the basis are chosen 
one at a time (and the corresponding vector to leave is determined by (4)). The basis transformation 
is then performed before the next vector to enter the basis is again chosen randomly. Since there are 
(n — m) nonbasis vectors, the random selection of a vector to enter the basis is performed by generating 
a random number from a uniform random distribution on the interval (0, 1), multiplying it by (n — m), 
and rounding the result to the next highest integer, say k. Then the A:th nonbasic vector is selected to 
enter the basis. The process is repeated 8 times, and the resulting extreme point is the first randomly 
generated extreme point within a distance 8 of Xo to be used in Step 2. 

An approach very similar to the above is proposed by D. Denzler [5]; the basic difference is that 
Denzler's method does not restrict its random search for an extreme point within a certain distance 
of the current best solution. Many iterations are thus wasted. Intuitively, it would seem that extreme 
points close to one which has a known good value of z should also yield fairly good values of z, and so 
it would appear to be a better procedure to restrict the search for better solutions to extreme points 
fairly close to Xo rather than attempt to search randomly the entire set of extreme points. 

2.4 The Cooper-Drebes and Cooper-Olson Heuristics 

Cooper and Drebes [3] have developed two other adjacent extreme point heuristics, which differ 
from the above three methods mainly in that in the Cooper-Drebes methods the objective function is 
modified at certain stages in the algorithms; in particular, the linear costs cj,j= 1, . . ., n, are modi- 
fied according to 

(8) c j = c j + ^,x j >0 

x j 

Cj=Cj, xj = 0, 

where "*/' implies the current value of variable Xj. At other stages in their heuristics, Cooper and 
Drebes remove from the basis the vector with the largest fixed charge (of the current basic variables) 
and insert into the basis the vector with the smallest fixed charge (of the current nonbasic variables). 
The two methods differ only in their decision rules for employing the various techniques discussed 
above. 

A modification of the Cooper-Drebes heuristics has been developed by Cooper and Olson [4]. The 
modification consists of randomly generating an extreme point a distance not more than 8 (where 
distance is defined as in Section 2.3) from the current local minimum, and then the Cooper-Drebes 
heuristic is again employed in an attempt to improve the solution. The random extreme point generation 
is called a "perturbation." The parameter 8 is allowed to vary from about m/3 to m/2. If the search 
following the last previous perturbation yields a local minimum worse than the best obtained so far, 
8 is reduced; if the objective function value is the same as the best, 8 is increased; if a better value is 
obtained, 8 is not changed. 



222 D. I. STEINBERG 

The computational experience reported on these methods [3, 4] indicates that they will yield the 
optimal solution a high percentage of the time and, moreover, can be expected to provide at least a 
near-optimal solution. 

2.5 A Separable Programming Approximation Technique 

Another technique which yields a local minimum is to formulate the fixed charge problem as a 
separable programming problem [7]. 

Letting Uj represent known upper bound values of Xj, consider the following representation of 
the objective function (1): 

(9) Z = J£ [djXy +(dj + CjUj)k 2j ] , 

where 

d°) xj = \ 2J Uj 

and 

(11) \0j+klj+\2j=l 

and where no more than two \kj can be positive for a given j, and these can only be positive if they are 
adjacent. The complete formulation becomes: 

N 

Minimize z= ]£ [djkij+ (dj + CjUj)\ 2 j], 
subject to 



^ aijUjkij — bu i = 1 , 



m 



AfcjSsO, k = 0, 1, 2 and j=l, . . ., n. 

When this problem is solved using a restricted basis entry Simplex method [7], the optimal solution 
will be a local minimum for the fixed charge problem. There is no evidence, however, that this procedure 
is any better than solving the fixed charge problem by using linear programming and ignoring the fixed 
charges. Moreover, no attempts are made to improve the local minimum obtained. 

3. A BRANCH AND BOUND ALGORITHM 

In attempting to develop an exact solution method for the fixed charge problem, it would be natural 
to begin by investigating the possibility of developing an adjacent extreme point method. However, 
since the objective function is not always monotonic (that is, there may not always exist a monotonically 
decreasing adjacent extreme point path to the optimal solution) some other approach is called for. In 
this section a branch and bound approach is presented. (A general knowledge of the branch and bound 
philosophy will be assumed in the following discussion. The reader is referred to [1] for a general dis- 
cussion of branch and bound methods.) 

The algorithm finds the global minimum solution without having to examine all basic feasible 
solutions of (2). In the next two sections, the general approach will be presented and the computational 
details will be given in Section 3.3. 



FIXED CHARGE PROBLEM 



223 



3.1 The Tree of Solutions 

Suppose the following tree is generated: 




The nodes are numbered from left to right and from top to bottom, so that each horizontal "line" 
of nodes may be thought of as a level; that is, the first level consists of only node 1, the second level 
consists of nodes 2 and 3, the third level consists of nodes 4, 5, 6, and 7, etc. The constraints (2) will 
be required to be satisfied at every node. In addition, each node (except node 1) will be required to 
satisfy some additional constraints, of the form "xj=0" or "xj > 0." Since only basic feasible (extreme 
point) solutions are of interest, the latter constraints actually require variable Xj to appear in the basic 
solution, at a positive level. At each node (except node 1) one additional constraint is imposed, so that 
the solution obtained at any given node k must satisfy this new constraint plus the constraints imposed 
at every other node in the path from node 1 to node k. Thus, a node at level L must satisfy {L— 1) 
additional constraints. 

If a path is terminated when the additional constraints imposed by the node are such that: 

1. There are m constraints of the form "xj > 0", or 

2. There are (n — m) constraints of the form "xj = 0", or 

3. No feasible solution exists when the current additional constraints are imposed, 

then it is clear that the collection of terminal nodes will contain all the basic feasible solutions to (2). 
That is, for terminal case 1, m ay's are required to be positive and hence a basis is uniquely determined. 
For case 2, (n — m) x/s must be zero, and hence the remaining xfs uniquely determine a basis; however, 
this basis may not yield a feasible solution, in which case terminal case 3 applies. 

Now, of course, the value of the objective function z can be computed for each extreme point and 
the minimal value obtained; however, it is desirable to reduce as much as possible the number of 
extreme points which must be examined. 



3.2 The Algorithm 

The following notations will be introduced: 



S* — {j\ x j must be positive at node k}, 
Sj = {j\xj must be zero at node k}, 
' z k — a lower bound of 2 at node k, and 
z* = the global optimal value of z. 



224 D - ! STEINBERG 

The constraint set at node k thus becomes: 

/ Ax = b 



(13) T k = 



x ^ 

Xj > 0, >eSf 

Xj = 0, jeS k 

x is an extreme point. 

Suppose that an upper bound estimate of z*, z, is obtained. Then a path can also be terminated 
at node k if 

(14) z k =s z, 

since the basic feasible solution corresponding to z has as good a solution as any of the solutions which 
could be obtained by continuing on the current path beyond node k. 

This procedure will cut down on the number of basic feasible solutions that need to be examined; 
its efficiency depends heavily on how good the bounds z k are and on how close z is to the true optimal 
value of z. 

Any heuristic method can be used to obtain an initial z\ thereafter, every time a path is terminated 
because a unique basic feasible solution has been obtained, the corresponding value of z, z, can be 
calculated and compared with the current z; if z is smaller than z, then z and its corresponding x 
replace z and its corresponding x as the current upper bound estimate of the optimal solution. A method 
for calculating z k is given in the next section. 

The branch and bound algorithm can be summarized as follows: 

1. An initial upper bound, z, of the optimal value of the objective function is obtained by means 
of a heuristic. This bound and its corresponding solution vector x are stored. 

2. The paths are then generated so that in branching from one node to another, the node at which 
the additional constraint is of the form "xj > 0" is selected as the next node in the path, until one of 
five path terminating rules occurs. The first three such rules are given on page 223. Rules 1 and 2 in- 
dicate that a unique solution has been reached (although it may not be feasible). Rule 3 pertains to 
infeasibility; that is, no basic feasible solution exists for the constraints (13). Rule 4 is the initial bound- 
ing procedure: A path is terminated if the lower bound z k of the objective function is as large as the 
current upper bound of z*; that is, if (14) holds, then the basic feasible solution corresponding to z 
has as good a solution as any of the solutions which could be obtained by continuing on the current 
path. Hence, the path is terminated. Rule 5 pertains to the fixed charges themselves. When the only 
fixed charges whose subscripts are not in either Si or S k are all zero, then as will be seen in the next 
section, the lower bound z k is the actual minimum value of the objective function (1) subject to the 
constraints (13). Hence, the path can be terminated (and if z k < z, then z k and its corresponding solution 
vector become the new upper bound). 

3. When the last path has been terminated, the current lower bound z and its corresponding x will 
be the optimal solution. 

3.3 Computation of the Lower Bounds 

In order to obtain the lower bound z k some initial computations are made. In particular, it will 
be useful to know the level of degeneracy that exists in (2); that is, how many basic variables can 
appear in a basic solution at a zero level. To obtain this information, the following n linear programming 
problems are solved: 



FIXED CHARGE PROBLEM . 225 

Minimize Zj = Xj, 

subject to: x is an extreme point of the convex set (2) and Xj is in the basis. 

The method for solving these restricted linear programming problems will be presented later in this 
section. 

If the optimal solutions to (15) are denoted by xf, respectively, then if all xf > 0, no degeneracy 
exists; if any x*—0, then the number of such variables determines the maximum level of degeneracy 
which could occur. 

It shall be assumed that: 

(16) di^d-2^ . . . ^d„. 

In practice (16) is really no restriction, since the vectors may only have to be reordered so that (16) 
holds. 

Let S$ = {j\j4S* US*}, 

P—{j\xj can be in the basis at a zero level}, 

c?=pn (Sfusf), 

N k = the number of elements in S k : , and 
N k i = the number of elements in Q. 

Thus S k is the set of subscripts whose corresponding variables do not appear as an additional 
constraint at node k; the set Q consists of those subscripts whose corresponding variables could appear 
in the basis at a zero level at node k. 

Consider the following problem: 

(17) Minimize z\~cx 

subject to: xeT k . 

The desired optimal solution to problem (17), for any node k, lies at an extreme point of the constraints 
(2) and also satisfies the additional constraints imposed at node k. Denote this optimal solution by 
z k and let 

(18) 2§ = 2 dj. 

Further, let z| denote the value of the sum of the last (m — N k — N k i ) fixed charges. Define 

(19) z k = z k + z k + z k . 

In order for z k to be a legitimate lower bound for z*, it is sufficient to show that for any solution xeT k 
(i.e., any solution x which satisfies the constraints at node k), 

(20) cx + dy^z k . 

But, 

cx + dy=cx+ V g?j + Y djy-j, 

(21) ;eSf Ui 

since for./€Sf, yj= 1 and for^'eS^', yj— 0. By definition, ex 3 s z k . 
(22) 



226 D. I. STEINBERG 

Since the maximum level of degeneracy possible at node k is N£, there must be at least (m — N£) fixed 
charges whose yj = 1 in any basic solution at node k; but,N k of these fixed charges are determined by S^ 
(and are included in z k ). Thus, there are at least (m — N k — Nq) remaining fixed charges whose yj must 
be unity. Since z k is defined to be the sum of the (tii — N* — Nq) smallest remaining fixed charges, it 
must be true that 

(23) 2 d M * 2 " 

its* 

Addition of the three relations (18), (22), and (23) yields the desired result (20). Hence z k as defined by 
Eq. (19) is a legitimate lower bound for the optimal solution which could be attained at node k. 

In order to complete the description of the algorithm it is only necessary to describe the procedure 
for solving the linear programming problems (15) and (17). Since (15) is merely a special case of (17) 
only the latter is considered below. 

Consider the following problem: 

(24) Minimize Z\ = ex 

(25) subject to: xeT k . 

The set of extreme points which satisfy the constraints (25) is a connected set; that is, an adjacent ex- 
treme point path exists between any two points in the set. Suppose an initial basic feasible solution to (25) 
is known, and that a Simplex algorithm is used to improve the objective function. Then, in a finite number 
of steps, a stage will be reached (call it Stage One) in which the only vectors which could be entered into 
the basis to yield a decrease in the objective function would require one of the vectors whose corresponding 
jeS k to leave the basis. (It is assumed throughout the discussion that none of the vectors whose corre- 
sponding jeS k are considered for entry into the basis.) Let the current solution at this stage be denoted 
by x, with corresponding i\ = ex. It is possible that some other extreme point £o exists such that zo < i\ , 
but that x and xo are not adjacent extreme points. In order to determine whether or not 5c is indeed an 
optimal solution to problem (17), the following constraint is appended to (25): 

(26) cx + x n + l = z l . 

The basis size is now increased to m + l, and the (m+ l)st basic variable for the current solution is 
x n + i = 0. Now, however, it is possible that one or more variables may be brought into the basis without 
removing any of the variables whose corresponding jeS k . If $ is not an optimal solution to (17), then a 
solution exists for which x n +i > 0. In order to determine whether such a solution does in fact exist, 
the following test for optimality is employed: 

Let R k denote the set of variables which are not in the basis for x (and whose j4S k ) . Then, for each 
xi in R k , first insert xi into the basis at a zero level, replacing x n + u and then enumerate the extreme 
points (of basis size m + 1 ) of the convex set determined by (25) and (26), by solving 

Maximize z = x,, 
subject to: Ax = b 

CX + Xn+l = Z\ 

(27) Xi — x„ +2 = xt 

x ^ 






FIXED CHARGE PROBLEM • 227 

xj > 0, jeS$ 
Xj = 0, jeS k 
x is an extreme point, 

where x, is the current value of x,-, initially zero. At each step in solving (27) the minimum increase in 
Xi is taken. In this manner every extreme point of (25) and (27) which contains x,- will be examined. If 
at any stage x n +i > 0, the procedure can be terminated, and by beginning with this new extreme point, 
the objective function (24) can be reduced in the usual way until Stage One is again reached. If, on the 
other hand, the optimal solution to (27) for each x\ is zero, or no extreme points for which x n + \ > have 
been found after solving (27) for all XitR k , then x is the optimal solution to problem (17). 

This procedure will not, in general, be as lengthy as it might at first appear. If x is the unique 
optimal solution to problem (17), then the problems (27) will all have optimal values of zero, and one 
iteration per problem will indicate this. Moreover, if x is not in fact optimal, is should not require very 
many iterations to find a solution for which x n +i > 0. 

It should also be noted that the true minimal solution to problem (17) is required only when 
i\ > (z — z k — z k ), where z is the upper bound estimate, since if z.\ =£ z the path will not be terminated 
and decreasing z\ further is of no interest. Hence the above test for optimality for problem (17) need 
be applied at any node only when i\ > z. 

The question still remains: How is an initial basic feasible solution to (17) obtained? The procedure 
is as follows: 

Let W\ = {xj\jeS k , x } not in current basis} 
W<i = {xj\J€S%, Xj > in current basis} . 

Thus, those XjeWi satisfy the constraints imposed at node k, and those XjeWt do not; in order to remove 
those XjtWt from the basis, the following linear programming problems may be solved, for each xfiW 2 : 

Minimize z=Xj, 
(28) subject to: Ax= b 

x^O. 

When z = 0, either xj has been removed from the basis, or it is in the basis at a zero level and can be 
removed; if in an optimal solution z > 0, no feasible solution exists with Xj not in the basis. 

It should be noted that once a variable XjeW 2 has been removed from the basis, then XjeW\\ from 
the theory of linear programming it is known that if it is desired that some variables not appear in a 
basic solution, they may merely be excluded from consideration for entry into the basis; thus, pro- 
ceeding to remove the variables XjeW 2 from the basis while never considering the variables XjeW\ for 
entry into the basis will yield in a finite number of steps a basic solution in which W\ = S^ and W 2 is 
empty. 

Now, the variables XjtS'l must be inserted into the basis: 

Let Fi = {xj\jeS\ and xj > 0} 
Y 2 = { Xj \J€S k 1 and Xj = 0}. 



228 



D. I. STEINBERG 



For those variables XjeY 2 the following linear programming problems are solved: 

Maximize z — Xj 

subject to: Ax = b 

(29) x ^ 

XjeYi remain in the basis 
at a positive level 

x is an extreme point. 

Since the XjeY x are already in the basis at a positive level, the problems (29) can be solved by the 
method discussed earlier in this section. If the optimal solution to any of the above problems is zero, 
then no feasible solution to problem (29) exists with that variable in the basis at a positive level; if any 
solution exists for which z > 0, then variable xj has successfully entered the basis (it is not necessary 
actually to optimize the objective function). At this point, Xj will now be in Yi and the next variable 
XjeY 2 can be inserted in the basis via problem (29). In a finite number of steps, a stage will be reached 
at which Y\ = S\ and Y 2 is empty. If none of the variables XjeS% were considered for entry into the basis 
during the procedure of solving the problems (29), then an initial basic feasible solution to problem 
(17) has been obtained. 

It should be noted that at no stage of the branch and bound procedure will more than one variable 
not satisfy the additional constraints imposed by any node k. Thus, the amount of actual computation 
necessary to obtain the initial solution is much less than one might at first suspect. 

Moreover, the amount of computer storage required by the algorithm does not build up as the 
number of nodes examined increases, but remains constant throughout the procedure (for a problem 
of any given size). This advantage does not exist in most other branch and bound algorithms. 

4. COMPUTATIONAL RESULTS 

The primary difficulty in finding fairly large problems to test heuristic methods is that the optimal 
solution must be known. Cooper and Drebes [3] have randomly generated 290 (5 X 10) problems with 

the following properties: 

(1) |ay|*£20, 

(2) l=£c jS £20, 

(3) 1 *£ dj *£ 999, 

(4) The density of A is 50 percent. 

The optimal solutions to these problems were obtained by total enumeration. 

In order to take advantage of these problems, larger problems were constructed by forming the 
A matrix as follows: 

""/I, • • • 

A 2 ■ • • 



A = 







A k 



FIXED CHARGE PROBLEM 229 

where the Ai are (5 X 10) matrices. The resulting problem is of size (5k X 10k) and its corresponding 
solution is known. Hence problems of varying size can be formulated to test the heuristic methods. 
Because of roundoff difficulties, some of the test problems are not included in the results. These are 
noted under the proper table. 

Since the heuristics tested are all adjacent extreme point methods, each iteration corresponds to a 
regular linear programming transformation. Thus, the primary means of comparison among different 
heuristic approaches is the number of iterations before the best obtained solution is found. Since these 
methods are not exact, the "best obtained solution" will not always be the optimal solution, and in 
some cases will be a function of the parameters related to the particular method. 

The approximate computer times (for an IBM 360/50) per iteration for the four groups of problems 
for both Heuristic One and Heuristic Two were (in seconds): 0.10, 0.18, 0.60, and 2.4. 

The Cooper and Drebes paper [3] reports computational experience with two heuristics: both 
heuristics were tested on a set of 253-(5 X 10) problems; 240 were solved exactly by the first method 
in an average of 75 total iterations per problem, and 245 problems were solved exactly by the second 
method in an average of 100 total iterations per problem. The three (5 X 10) problems which neither 
heuristic solved exactly were all less than 10 percent from optimality. No information on the number 
of iterations required to achieve the best solution is given. 

On a set of 70-(15 X 30) problems, the two methods were combined; the resulting heuristic solved 
63 problems exactly with an average of 1200 iterations per problem. The nonoptimal problems were 
less than 3 percent from optimality. 

4.1 Heuristic One 

Problems of size (5 X 10) ,( 15 X 30) , (30 X 60) , and (50 X 100) were used in the testing of Heu- 
ristic One. Of the 250-(5 X 10) problems, 235 were solved exactly with olq = 2 and /3o = 2. The average 
number of iterations required to achieve the final solution was 3 and the average total number of 
iterations per problem was 6.8. Table 1 contains a summary of the 15 nonoptimal problems. When ao 
was increased to 3 and /3o to 4, problems 2, 94, and 128 were solved exactly. No further increase in 
ao and /3o yielded any improvements in the 12 remaining suboptimal problems. 

Ninety (15 X 30) problems were attempted, and with a = 2, /3 — 2, the optimal solutions to 75 
problems were obtained, with an average of about 5 iterations before the best solution was obtained 
and an average of 12 total iterations per problem. No further improvement resulted when ao was 
increased to 10; the average total number of iterations increased to 16. /3o was increased again to 15, 
with an improvement in only one of the suboptimal solutions-problem 43. The optimal solution is 
Z* = 8317.7 and the best previously obtained solution was Z = 9172.4; with /3 = 15, the solutions 
8964.7, 8593.0, and 8525.8 were obtained. The average total number of iterations was 23. Table 2 
contains a summary of the suboptimal (15 X 30) problems. 






230 



D. I. STEINBERG 



Table 1. Heuristic One-(5 X 70) Suboptimal Problems 



Problem 


Optimal 


Best obtained 


Percent 


number" 


solution 


solution 


from optimal 


2 


2830.4 


"3317.6 


14.7 


57 


4015.8 


4017.9 


0.05 


70 


3363.0 


3468.3 


3.1 


86 


2359.4 


2380.7 


0.9 


93 


2351.7 


2436.6 


3.6 


94 


2373.0 


"2899.1 


1.1 


99 


2418.0 


2473.0 


2.3 


128 


3459.0 


» 3470.6 


0.3 


131 


1543.0 


1766.1 


14.5 


137 


3593.1 


3885.1 


8.1 


181 


2826.2 


3183.0 


12.6 


188 


2417.1 


2440.1 


0.95 


235 


2080.9 


2155.9 


3.6 


267 


3072.2 


3170.9 


3.2 


279 


2323.4 


2533.7 


9.1 



"Problems not attempted: 29, 36, 51, 55, 66, 67, 73, 75, 82, 88, 98, 117, 119, 120, 126, 130, 134, 141, 147, 158, 162, 164, 166, 
176, 185, 201, 204, 206, 211, 215, 219, 226, 236, 237, 238, 248, 250, 260, 264, 266, 275. 
b Exact solutions obtained with a = 3 and /3o = 4. 



Table 2. Heuristic One-(15X30) Suboptimal Problems 



Problem 


Optimal 


Best obtained 


Percent 


number" 


solution 


solution 


from optimal 


20 


7676.4 


7726.9 


0.67 


29 


7522.6 


7544.0 


0.28 


30 


9452.9 


9458.2 


0.06 


31 


9201.7 


9286.5 


0.89 


32 


6666.4 


7192.6 


7.89 


33 


8738.3 


8793.3 


0.63 


40 


11,117.3 


11,365.9 


2.23 


43 


8317.7 


8525.8 


2.50 


45 


6785.6 


6938.3 


2.25 


46 


8517.7 


8809.7 


3.43 


53 


6497.6 


6647.7 


2.31 


61 


7651.0 


8007.7 


4.66 


63 


7836.0 


7859.6 


0.29 


85 


6307.2 


6635.2 


5.20 


93 


8592.7 


9047.5 


5.29 



" Problems not attempted: 1, 44, 64, 68, 79, 89. 



FIXED CHARGE PROBLEM 



231 



Thirty-six (30x60) problems were attempted and 26 were solved optimally with a — 6 and /3 = 25. 
The average number of iterations before the best solution was obtained was about 11 iterations and the 
total number of iterations per problem averaged about 42. A summary of the suboptimal (30X60) prob- 
lems is given in Table 3. 

TABLE 3. Heuristic One-(30 X 60) Suboptimal Problems 



Problem 


Optimal 


Best obtained 


Percent 


number " 


solution 


solution 


from optimal 


15 


16,975 


17,002 


0.04 


16 


15,868 


16,479 


3.98 


17 


18,441 


18,4% 


0.30 


23 


15,303 


15,748 


2.91 


25 


16,886 


16,987 


0.50 


27 


13,960 


14,110 


1.07 


31 


14,944 


15,301 


2.39 


32 


17,350 


17,440 


0.52 


37 


14,716 


15,181 


3.16 


43 


13,946 


14,274 


2.35 



° Problems not attempted: 1, 5, 10, 20, 22, 24, 34, 40, 45. 
Nine-(50X100) problems were also tested, with a ~ 10 and /3o = 40; the results are given in Table 4. 



Table 4. Heuristic One-(50 X 100) Problems 



Problem 
number 


Optimal 
solution 


Best obtained 
solution 


Percent 
from optimal 


No. of iterations 

to best obtained 

solution 


Total No. of 
iterations 


1 
2 
5 
12 
14 
15 
18 
19 
20 


26,761 
28,842 
25,545 
29,350 
24,444 
29,380 
26,844 
25,903 
27,822 


26,761 
28,842 
25,545 
30,182 
25,085 
29,481 
26,844 
26,350 
27,899 




47 
22 
21 
14 
28 
13 
18 
24 
18 


87 
62 
61 
49 
46 
48 
51 
49 
48 






2.83 
2.57 
0.34 


1.73 
0.28 


Average No. of Iterations 


22.8 


55.7 



4.2 Heuristic Two 

The same sets of problems used to test Heuristic One were employed in the testing of Heuristic 
Two. In some cases, again, roundoff difficulties caused some problems to be discarded, and these 
problems are again noted in the tables. 



232 



D. I. STEINBERG 



Of the 268- (5X10) problems attempted by Heuristic Two, 255 were solved optimally. A summary 
of the suboptimal problems appears in Table 5. The average number of iterations required to achieve 
the final solution was 6, and the average total number of iterations was 15.2. 

TABLE 5. Heuristic Two-(5 X 10) Suboptimal Problems 



Problem 


Optimal 


Best obtained 


Percent 


number " 


solution 


solition 


from optimal 


28 


2993.4 


3025.4 


1.1 


50 


2383.9 


2813.4 


18.0 


71 


3058.3 


3470.5 


13.5 


93 


2351.7 


2436.6 


3.6 


99 


2418.0 


2473.0 


2.3 


108 


2414.8 


2423.8 


0.4 


119 


3467.6 


3560.6 


2.7 


120 


3579.1 


3734.5 


4.3 


131 


1543.0 


1766.2 


14.5 


137 


3593.1 


3885.1 


8.1 


154 


2589.5 


2622.7 


1.3 


245 


3574.7 


3635.8 


1.7 


279 


2323.4 


2533.7 


9.1 



"Problems not attempted: 23, 27, 30, 52, 61, 67, 76, 81, 94, 98, 110, 128, 149, 158, 165, 202, 210, 224, 225, 239, 249, 254. 

Of the 287- (5X10) problems attempted by at least one of Heuristic One and Heuristic Two, 280 
were solved optimally; the following 7 problems were not solved optimally by either method: 93, 99, 
119, 120, 131, 137, and 279. 

Eighty-four (15X30) problems were attempted by Heuristic Two, and 74 were solved optimally. 
The average number of iterations required to obtain the final solution was 12, and the average total 
number of iterations per problem was 43. 

Table 6. Heuristic Two-(15 X 30) Suboptimal Problems 



Problem 


Optimal 


Best obtained 


Percent 


number ° 


solution 


solution 


from optimal 


1 


8096.2 


8142.6 


0.57 


31 


9201.7 


9286.5 


0.89 


36 


6162.5 


6171.6 


0.15 


40 


11,117.5 


11,365.9 


2.23 


43 


8525.8 


8537.4 


2.50 


44 


6807.3 


6905.5 


1.44 


46 


8517.7 


8809.7 


3.43 


52 


7880.0 


8214.1 


2.97 


79 


7765.0 


7840.0 


2.45 


93 


8592.7 


8803.0 


5.29 



"Problems not attempted: 16, 18, 21, 22, 26, 27, 33, 37, 61, 72, 80, 83. 



Of the 96 (15 X 30) problems attempted by at least one of Heuristic One and Heuristic Two, only 
the following 5 problems were not solved optimally: 31, 40, 43, 46, and 93. 



,u i^riAKGE PROBLEM 



233 



TABLE 7. Heuristic Two-(30 X 60) Suboptimal Problems 



Problem 


Optimal 


Best obtained 


Percent 


number 


solution 


solution 


from optimal 


1 


15,490 


15,536 


0.30 


20 


18,407 


18,655 


1.34 


22 


15,333 


15,535 


1.32 


23 


15,303 


15,595 


1.91 


26 


16,848 


17,346 


2.91 


29 


15,773 


15,889 


0.73 


34 


16,912 


17,076 


0.97 


41 


15,601 


15,662 


0.31 


45 


14,243 


14,267 


0.17 


47 


15,852 


16,062 


1.32 

i 



"Problems not attempted: 5, 9, 11, 13, 14, 16, 17, 19, 31, 36, 40, 42, 43. 

Thirty-six (30 X 60) problems were attempted by Heuristic Two, and 26 were solved optimally. 
The average number of iterations required to obtain the final solution was 25, and the average total 
number of iterations per problem was 81. 

Of the 46 (30 X 60) problems attempted by at least one of Heuristic One and Heuristic Two, 42 
were solved optimally; the following problems were not solved optimally by either method: 1, 23, 34, 
and 45. 

Heuristic Two was not tested on any larger problems. 



4.3 The Branch and Bound Method 

Fifteen (5 X 10) problems and 10 (15 X 30) problems were solved by the branch and bound pro- 
cedure. The results are given in Tables 8 and 9. The total possible number of nodes that might have to be 

Table 8. The Branch and Bound Algorithm-(5X 10) Problems 



Problem 


Number 


Time 


number 


of nodes 


(minutes) 


1 


66 


0.23 


2 


24 


0.16 


3 


20 


0.16 


4 


44 


0.26 


5 


20 


0.13 


6 


14 


0.11 


7 


20 


0.13 


8 


74 


0.31 


9 


10 


0.09 


10 


38 


0.23 


11 


34 


0.17 


12 


12 


0.10 


13 


40 


0.16 


• 14 


24 


0.13 


15 


42 


0.30 


Average 


32 


0.16 



234 



D. I. STEINBERG 
TABLE 9. The Branch and Bound Algorithm-(l 5X30) Problems 



Problem 


Number 


Time 


number 


of nodes 


(minutes) 


6 


1242 


17.9 


8 


672 


9.1 


9 


3144 


47.3 


14 


1134 


18.3 


18 


308 


8.5 


29 


2294 


43.0 


34 


502 


26.4 


35 


334 


4.8 


37 


904 


11.6 


38 


1546 


24.4 


Average 


1208 


21.1 



examined is 2" +1 — 1. The average number of nodes that were examined was 32 and 1208, for the two 
groups of problems, respectively. Also included in the tables is the computer time required to move 
through the nodes; these times are for an IBM 7072. It is felt that these times could be greatly improved 
if the computer code were made more efficient. There are many improvements that could be made in 
the code, such as rewriting it in machine language and utilizing a more efficient Simplex method (such 
as revised Simplex). 

The average number of nodes examined for the (5 X 10) problems was 2 5 out of a possible total 
of approximately 2 10 nodes; for the (15 X 30) problems the average number of nodes actually examined 
was approximately 2 10 out of a possible total of approximately 2 30 . Hence, the ratio of the exponents 
of the number of nodes examined to the number of total possible nodes decreased when the number 
of constraints was tripled (from 5 to 15). If it is assumed that this ratio continues to decrease linearly, 
then the number of nodes that would have to be examined for a problem of size (45 X 90) is about 
2 20 nodes. 

Hence, it appears that the number of nodes examined grows rapidly with the number of variables, 
and so the efficiency of the computer code is extremely important. 

5. SUMMARY AND CONCLUSIONS 

The computational results of Section 4 clearly indicate that the heuristic algorithms are orders 
of magnitude faster than the branch and bound algorithm. Moreover, the heuristics almost always 
provide very good, if not optimal, solutions. When two or more heuristics are combined, the likelihood 
of finding the optimal solution increases considerably. In fact, when the two heuristics of Cooper and 
Drebes [3] are combined with Heuristic One and Heuristic Two, all 290 of the (5 X 10) problems were 
solved optimally, only one problem of the 96 (15 X 30) problems (problem 53) was not optimally solved. 

The branch and bound algorithm presented herein is a new, exact algorithm for the solution of 
the fixed charge problem and appears to be much better than the mixed integer-continuous variable 
approach discussed in Section 1.2. Such an exact method of solution is essential for those situations 
in which the true optimal solution must be determined. 

It is felt that many improvements could be made in the branch and bound algorithm, in addition 
to refinements in the computer code used to test the algorithm. It is hoped that a better method of 



FIXED CHARGE PROR.EM 235 

determining the true maximum level of degeneracy of the convex set of solutions can be found, so that 
a better lower bound for the objective function at any given node may be computed. Other refinements 
in the handling of the degeneracy problem could probably also be found. 

It is suggested that research be conducted on the problem of generalizing the branch and bound 
approach described herein, so that it would be applicable to more general separable concave functions. 
Another generalization which should be considered is the case in which the fixed charges are dependent 
upon more than one variable. For example, a fixed charge dtj could be incurred only if the corresponding 
variables xt and Xj are both positive. 

6. ACKNOWLEDGMENTS 

The author wishes to acknowledge the support of the Washington University Computing Facilities 
through NSF Grant G-22296 and also the Atomic Energy Commission for partial support of the research 
under Research Contract No. AT(11-1)-1493. 

REFERENCES 

[1] Agin, Norman, "Optimum Seeking with Branch and Bound," Management Science, 13 B- 176- 185 (Dec. 1966). 

[2] Balinski, M. L., "Fixed-Cost Transportation Problems," Nav. Res. Log. Quart. 8 41-54 (Jan. 1961). 

[3] Cooper, L. and C. Drebes, "An Approximate Solution Method for the Fixed Charge Problem," Nav. Res. Log. Quart. 14, 

101-113 (Mar. 1967). 
[4] Cooper, L. and A. M. Olson, "Random Perturbations and the MI-MII Heuristics for the Fixed Charge Problem," Department 

of Applied Mathematics and Computer Science, Washington University, Report No. COO-1493-7. 
[5] Denzler, D. R., "An Approximative Algorithm for the Fixed Charge Problem," Nav. Res. Log. Quart. 16, 411-416 (Sept 1969). 
[6] Hadley, G., Linear Programming (Addison- Wesley Publishing Company, Inc., Reading, Massachusetts, 1962). 
[7] Hadley, G., Nonlinear and Dynamic Programming (Addison- Wesley Publishing Company, Inc., Reading, Massachusetts, 

1964). 
[8] Robers, P. and L. Cooper, "A Study of the Fixed Charge Transportation Problem," Department of Applied Mathematics and 

Computer Science, Washington University, Report No. COO- 1493-9. 






ON COMPUTING THE EXPECTED DISCOUNTED RETURN IN A 

MARKOV CHAIN* 



David F. Hitchcock 

and 
James B. MacQueen 






University of California, Los Angeles 



ABSTRACT 

The discounted return associated with a finite state Markov chain X t , X 2 ■ ■ ■ is given by 
g(Xi) + ag(X 2 ) + a 2 g(X 3 ) + . . ., where g(x) represents the immediate return from state 
x. Knowing the transition matrix of the chain, it is desired to compute the expected discounted 
return (present worth) given the initial state. This type of problem arises in inventory theory, 
dynamic programming, and elsewhere. 

Usually the solution is approximated by solving the system of linear equations character- 
izing the expected return. These equations can be solved by a variety of well-known methods. 
This paper describes yet another method, which is a slight modification of the classical itera- 
tive scheme. The method gives sequences of upper and lower bounds which converge mono- 
tonely to the solution. Hence, the method is relatively free of error control problems. 

Computational experiments were conducted which suggest that for problems with a large 
number of states, the method is quite efficient. The amount of computation required to obtain 
the solution increases much slower with an increase in the number of states, N, than with 
the conventional methods. In fact, computational time is more nearly proportional to /V 2 , 
than to N :i . 

I. Introduction 

Consider a Markov chain X x , X 2 , . . ., taking values in the finite set S= [1,2, . . ., /V] that has 
transition probabilities p(y;x) = Pr[X t+ i — y\Xt = x]. When in state x an amount g{x) is received 
immediately. The discounted return associated with passing through the sequence of states 
X\, X2, ■ ■ -, is 

g{X l ) + ag(X i ) + c?g(X 3 ) + . . ., 
where a is the discount factor (0^a<l). 

The expected discounted return given Xi = x, represented by u(x) (also called the present worth 
of state x), is well known to be the unique solution of the equation 

(!) u(x) = g(x) + a^ u(y)p(y;x). 

An efficient machine algorithm for solving (1) is desired. 

The problem of computing the expected discounted return frequently arises in inventory problems, 
replacement problems, and elsewhere (c.f. [3], [4]). 

If N (the number of states in S) is reasonably small, then (1) can be solved by familiar techniques, 
such as Gausian elimination or the simplex method, and is a straight forward problem in machine com- 
putation. However, if there is a large number of states, the amount of computation, which is generally 

*This work was supported by the Western Management Science Institute under a grant by the Ford Foundation, and by 
the Office of Naval Research under Contract No. 233(75), Task No. 047-003. Computational facilities were provided by the 
Western Data Processing Center, University of California, Los Angeles. 

237 



238 



D. F. HITCHCOCK AND J. B. MACQUEEN 



proportional to N 3 (cf.[l, p. 215]) using such methods, may be prohibitively expensive. Also, with these 
methods error control problems can arise even in problems of moderate size. 

This note describes some computational experiments using a slightly modified version of the 
classical iterative method for computing u. The method determines sequences of functions {u' n } and 
{u'n) on S which converge monotonely to u from below and above, respectively. Then u n = (u' n + u'^)/2 
is taken as the estimate of u. The experiments indicate that for "large" problems (say /V ^ 100) , the 
amount of computation required to obtain a solution whose percentage error at each state is less than 
a preassigned figure, increases much more slowly than /V 3 . For example, in a number of the experiments 
with TV = 500, approximately 16(500 2 ) elementary (multiply and add) computations were required to 
achieve 0.005-percent accuracy. As with the estimate produced by the classical iteration scheme 
(see Eq. (5) below) convergence is found to be exponential, but at a rate significantly faster than the 
classical method. The monotonely convergent upper and lower bounds, which are easily calculated, 
provide a very satisfactory means of controlling error due to round-off, which may accumulate with 
other methods. 

The method starts with any function v\ which is at a selected state seS, with the functions t>2, 
t>3, . . ., being defined by 

(2) v„ + 1 (x) = g(x) + a 2 Vn(y)p(y; x) - (g(s) + a J v n (y)p(y; *) )• 

Thus v n (s) = 0, n = 1, 2, . . 



(3) 



. ... Let 

L' n = min x 
L" n = max x 



g(x) + a ^ v n (y)p(y; x) - v„(x) 

y 
g(x) + a ^ v n {y)p(y; x) - v„(x) 



The sequences {u' n } and {u" n } are then given by 

(4) u' n = v„+ (1- a)~ l L' n 

u: = v n + (l-a)-'L; 

The sequences {u' n } and {u'^\ have the properties: 

(i) u' n ss u =£ ul, (ii) u' n s£ u' n + 1 , u" n ^ u'„ +u and (hi) u' n -* u, u" n -* u. 

These properties are obtained directly as special cases of Theorem 2 in [2]. 

The modification represented by the calculation of u„ is to be compared with the classical iterative 
method. The latter is defined by an initial function U\ and the rule 



(5) 



u„ + i(x) = g(x) + a^ u n (y)p(y; x). 



If ui{x) =vi(x) + «i(s), then v n (x) = u n {x) — u n (s), and 

«„U) = «„(x) + V2(i-a)-'(y;+y;'), 

y^ = min x g(x) + a^u n (y)p(y; x) — u n {x) 
L v 

Jl=m<ix x g{x) +a^u n (y)p(y; x) —u„(x) 



(6) 
where 

(7) 



COMPUTING RETURN IN MARKOV CHAIN 239 

II. Experimental Procedure 

To explore the computational efficiency of the suggested method, it was decided to investigate 
four main types of transition matrices, with TV equal to 10, 50, 100, 500, or 1,000. The discount factor 
was varied over the four values 0.90, 0.95, 0.99, and 0.995 representing the likely range of applications. 
In all the experiments, the initial function v\ was identically zero. 

The four main types of transition matrices were as follows: 

(I) Random. In this type of matrix, the row entries were randomly generated by selecting for each 
entry a uniformly distributed random number between and 1 and then normalizing by dividing by 
the row total. 

(II) Large probability in the left hand column. Here the first column has a specified value po (i.e., 
p(l; x) —po for x= 1,2, . . ., N). The value of po was varied over the three values 0.1, 0.5, and 0.9, 
with the remaining entries generated randomly as above (summing to 1 — po, of course). 

(III) Large probability on the diagonal. Here all diagonal elements were set equal to po, and the other 
entries were random as before. The value of po was varied as above. 

(IV) Large probabilities scattered. In this type of matrix, the first row was generated by taking four 
probabilities pi=p /2, p 2 = po/4, p 3 = po/8, and p 4 = p /8 and spacing them evenly in the first row. 
For example, for /V=100 we have p(l; l)=pi, p(25;l)=p 2 , p(50; l)=p 3 , and p(75;l)=p 4 . The 
remaining rows of p(y; x) were obtained from the first row by cyclic permutations, the last entry in 
the first row becoming the first entry in the second row, etc. 

It was decided arbitrarily to have g(x) vary in multiples of 5 as follows: for n—10, g(l) = 5, 
g{2) = 10, . . ., g(10) = 50. For the remaining values of N, all of which were multiples 10, the first 
N/10 values of g were 5, the second N/10 values were 10, etc., with the last A710 values being 50. Evi- 
dence that this method of choosing g, while arbitrary, was not malicious, was obtained by repeating 
all the experiments for N= 100 for matrices of Type III (large probabilities on the diagonal) with a 
different scheme for choosing g. Each value of g(x) , x— 1, 2, . . ., N, was determined as a uniformly 
distributed random number between and 100. Matrices of Type HI were chosen for this additional 
experimentation since they were the most difficult of the four in terms of convergence times. The 
number of iterations (see below, Table 1) using the randomly determined g was essentially identical 
to those obtained with the fixed g, differing by 0, 1 or at most 2 iterations. 

III. Results 

The main results are presented in Table 1. The entries in the table represent the number of 
iterations (the value of n in Eq. (4)) required for the maximum possible percentage error (MPPE) of 
the estimate u n to be less than 0.005 or 0.5 percent. Thus, for the indicated value of n, max x \u(x) 
— iin(x) | /u(x) is less than 0.00005 or 0.005. The upper and lower entries in each cell in the table 
correspond respectively to the 0.005 and 0.5 percent figures. The values of n were actually computed 
by allowing the computer program to run until max x {u'n{x) — u' n (x))lu' n {x) =£0.0001, which implies 
max x |«(x) — u(x)\lu(x) =£ 0.00005. The 0.5-percent values were then determined by inspection of 
the output which included the MPPE for each iteration. 

It should be noted that the MPPE may be considerably larger than the actual percentage error. 
As a typical example, consider the matrix of Type III with P = 0.90, /V= 100 and a = 0.90 where the 
MPPE after 16 iterations was 6.4 percent. However, ui 6 (x) differed by at most 1.4 percent from u 52 (jc), 
and the latter (which corresponds to the last iteration performed) has at least 0.005 percent accuracy. 



240 D - F HITCHCOCK AND J. B. MACQUEEN 

The MPPE would not be a meaningful criterion of accuracy if u(x) were very near zero for some x. 
This possibility was avoided by the choice of strictly positive values of g. The MPPE can also be 
rendered meaningless by adding a large positive constant to each value of g, since if u is the solution 
with immediate return function g, u+A/(l — a) is the solution with immediate return g + A. Thus by 
making A sufficiently large, the MPPE in the approximate and easily computed estimate, 0. — A/(l — a) 
can be made arbitrarily small. Again, the choice of g avoids this trivial extreme. 

Several general findings were noted. First, the MPPE exhibited a strong tendency to decrease 
exponentially with the number of iterations, i.e., as kr n ,0 > r > 1. Moreover k is typically much less 
than unity. This was clearly evident on plotting the error against n on semi-logarithmic paper. The 
maximum percentage error for the classical method defined by equation (5), also decreases exponen- 
tially, approximately as Ka" where K is near unity. However, r was appreciably less than a, the ratio 
r/a being .9 or less in all the experiments. Thus the ratio of percentage errors for the two methods after 
30 iterations, say, would be around .05 and perhaps smaller. 

A second finding is that the method automatically takes advantage of certain simplifying features 
of the transition matrices. Thus in matrices of Type I and II, the distribution of the next value of g 
is very similar from state to state, and the sequence g(Xi), g(X 2 ), . . ., is a sequence of similarly 
distributed, nearly independent random variables. For that reason, a good estimate for u(x) is simply 
the immediate return g( x) plus the mean value of g times a/ (I — a), and a good computational scheme 
ought to take advantage of this simplicity. The experimental results show that the method does. Some- 
what surprisingly, when the rows of the transition matrix are identical, then 

U2(x) =u'i(x) =g(x) +a(l-a)- 1 ^g(y)p(y; x)=u(x), 

y 

and the exact result is obtained immediately (given V\ — 0). 

This suggests a general rule, that the less the effect of the initial value of x on the sequence of 
ensuing returns, the less the amount of computation required. On this basis, one would predict the 
Type IV matrix to be easier to solve than the Type HI matrix, which in fact is the case. 

A third resul suggested by Table 1 is that the amount of computation required to achieve a given 
MPPE increases approximately as N 2 . There are approximately N 2 elementary computations (multiply 
and add being counted as an "elementary computation") per iteration, so that the amount of computa- 
tion required in n iterations is proportional to nN 2 . As is evident in Table 1, for a given combination of 
matrix type, Po and a, n increases slowly and at a decreasing rate, or does not increase at all, as /Vis 
increased. The worst case from the point of view of the relative increase in computation with N is the 
Type IV matrix with a = 0.95 and po = 0.9. However, even in this instance only nine more iterations 
are required to achieve 0.5-percent error as N is increased by a factor of 10 from 50 to 500. Note, how- 
ever, that in Table 1 with N — 10, there are a good number of instances where more than 10 iterations 
and hence more than N 3 elementary computations are required. Thus the relative advantage of the 
method as against conventional method requiring approximately /V 3 elementary computations, would 
appear to hold only for large values of N. 

Finally it is to be noted that the suggested method seems to be insensitive in an appropriate way 
to near zero entries in the transition matrices. This is perhaps a major reason for the rapid convergence. 
Type I matrices for which 0.005-percent accuracy was achieved in four iterations, particularly illustrate 
this result as it is obvious that for N = 500 or N = 1 ,000 the entries are extremely small. 



COMPUTING RETURN IN MARKOV CHAIN 241 

TABLE I. Iterations Required to Achieve 0.5 and 0.005 percent Maximum Possible Error 



Matrix Type 


I 


II 


III 


IV 


\ N 
a \ 


10 


50 


100 


500 


1000 


Prob 
Po 


\ N 
a \ 


10 


50 


100 


500 


1000 


\ N 
\ 


10 


50 


100 


500 


1000 


\ N 
a \ 


50 


100 


500 


1000 


0.995 


7 

4 


5 

4 


5 
3 


4 
3 


4 
3 


0.9 


0.995 


4 
3 


4 
3 


4 
3 


4 
3 




0.995 


83 
46 


90 

51 


91 
51 






0.995 


58 
30 


61 
32 








49 




32 










0.99 


7 
4 


5 
4 


5 
3 


4 
3 


4 
3 


0.99 


4 
3 


4 
3 


4 
3 


4 
3 




0.99 


80 
45 


87 
49 


87 
49 






0.99 


43 
19 


54 
24 


























0.95 


7 
4 


5 
4 


5 
3 


4 
3 


4 
3 


0.95 


4 
3 


4 
3 


4 
3 


4 
3 


3 
3 


0.95 


65 

37 


66 
38 


69 
39 






0.95 


36 
16 


43 

20 






37 




25 








0.90 


6 

4 


5 
4 


5 
3 


4 
3 


4 
3 


0.90 


4 
3 


4 
3 


4 
3 


3 
3 


3 
3 


0.90 


50 
28 


52 
30 


52 
29 






0.90 


25 
12 


37 
17 






29 




20 










0.5 


0.995 


5 
3 


5 
3 


5 
3 


6 
3 




0.995 


15 
9 


16 
9 


16 
9 






0.995 


15 
9 


19 
10 


























0.99 


5 
3 


5 
3 


5 
3 


5 
3 




0.99 


15 
9 


16 
9 


16 
9 






0.99 


11 
6 


12 

7 


























0.95 


5 
3 


5 
3 


5 
3 


5 
3 


4 
3 


0.95 


15 
8 


15 
9 


16 
9 






0.95 


11 
6 


12 

7 


12 

7 




9 










0.90 


5 
3 


5 
3 


5 
3 


4 

3 


4 
3 


0.90 


15 
8 


15 
9 


14 
9 






0.90 


10 
6 


11 
6 


12 

7 




8 










0.1 


0.995 


7 
4 


6 
3 


6 

3 


5 
3 




0.995 


7 
4 


6 

4 


6 
4 






0.995 


6 
4 


6 
3 


























0.99 


7 
4 


6 
3 


6 
3 


5 
3 




0.99 


7 
4 


6 
4 


6 

4 






0.99 


6 
4 


5 
3 


























0.95 


7 
4 


5 
3 


5 
3 


5 
3 


4 
3 


0.95 


7 
4 


6 

4 


6 
4 






0.95 


6 
4 


5 
3 


5 
3 




4 










0.90 


7 
4 


5 
3 


5 
3 


4 
3 


3 
3 


0.90 


7 
4 


6 

4 


6 

4 






0.90 


5 
3 


5 
3 


5 
3 


4 














4 


4 


3 



REFERENCES 



[1] Kunz, Kaiser S., Numerical Analysis (McGraw-Hill, New York, 1957). 

[2] MacQueen, James B., "A Modified Dynamic Programming Method for Markovian Decision Problems," J. Math. Anal, and 

Appl. 14, 38-43 (1966). 
[3] Scarf, Herbert E., "A Survey of Analytical Techniques in Inventory Theory." In Multistage Inventory Models and Techniques, 

Herbert E. Scarf, Dorothy M. Gilford, and Maynard W. Shelly, editors (Stanford Univ. Press, 1963). 
[4] Wagner, Harvey M., Michael O'Hagen, and Bertil Lundt, "An Emperical Study of Exactly and Approximately Optimal In- 
ventory Policies," Management Science, 11, 690-723 (1965). 



MORE ADO ABOUT ECONOMIC ORDER QUANTITIES (EOQ) 



Victor J. Presutti, Jr. 

and 

Richard C. Trepp 

Operations Analysis Office 

Headquarters, Air Force Logistics Command 



ABSTRACT 



This paper is concerned with the determination of explicit expressions for economic 
order quantities and reorder levels, such that the cost of ordering and holding inventory is 
minimized for specific backorder constraints. Holding costs are applied either to inventory 
position or on-hand inventory, and the backorder constraint is considered in terms of the total 
number of backorders per year or the average number of backorders at any point in time. 
Through the substitution of a new probability density function in place of the normal p.d.f., 
explicit expressions are determined for the economic order quantities and the reorder 
points. The resulting economic order quantities are independent of all backorder constraints. 
It is also concluded that under certain conditions, the minimization of ordering costs and 
inventory holding costs (applied to inventory position), subject to a backorder constraint, 
is equivalent in terms of reorder levels to minimization of the safety level dollar investment 
subject to the same backorder constraint. 

INTRODUCTION 

We are concerned with a single-echelon, multi-item, continuous review inventory system in which 
the process generating demands does not change with time. When an item's reorder level is reached, 
an order is placed. After the elapse of a leadtime, which is assumed constant, the entire order is re- 
ceived. For each item, it is desired to determine the constant reorder quantity and reorder point which 
will minimize (summing all items in the inventory) the holding and procurement costs subject to an 
overall backorder constraint. The holding cost may be a function of either on-hand inventory or inven- 
tory position. We will constrain total units backordered or average number of units in a backorder 
position, which, as is well known, is also the unit-years backordered per year. In either case, the 
backorders may be weighted by item essentiality. 

MODEL DERIVATION 
Definitions 

The following definitions apply throughout the discussion: 
Q: order quantity (units) 
r: reorder level (units) 
D: annual demands (units) 
k: safety factor 

cr: standard deviation of unit demand during a leadtime 
a: holding cost factor 
A: ordering cost (dollars per order) 
c: item cost (dollars per unit) 
fju: mean leadtime demands 
Z: item essentialtiy (relative military worth) 

243 



244 



V. J. PRESUTTI, Jr. AND R. C. TREPP 



x: units demanded during leadtime 

T)\ constraint — total essentiality-weighted units backordered 

/3: constraint — expected number of essentiality-weighted units in a backordered position 

at any point in time 
t: leadtime. 

We also require: 
i) r = fx + kcr 

ii)/(x),the probability density function (p.d.f.) describing the distribution of units demanded 
during the leadtime: 



(1) 



2^ ex P ( 



X — /JL 



for all x. 



It is easily shown that this implies: 



(2) 
and 

(3) 



Pr[x > fx+ ko-]= 0.5 exp (- V2k) for k 2* 



Pr[x > fJL + k<r] = 1 - 0.5 exp ( V2A) for k < 0. 



The form of f(x) was inspired by the work of L. L. Parker [3] who essentially assumed demands 
were normally distributed, and then approximated the complementary cumulative function by a function 
of the form ae~ bx . We originally chose a and b such that the mean absolute deviation from the normal 
distribution over the range =s k =s 3 was minimized (a = 0.57, 6=1.35). We now consider, however, 
that a = 0.5 and b— v2 is more appropriate.* In addition to defining a p.d.f. that is symmetric about 



the mean (/u,), these values also cause Pr 



hi < ?—£ < k 2 



to be independent of fi and u (a desirable 



attribute also shared by the normal p.d.f.). Further, with these values, the expected backorders at a 
random point in time for k = is the same as if demands were normally distributed, for large Q/cr. This 
similarity was considered desirable since we are primarily interested in average backorders (rather than 
total backorders) and we will often have no safety level. 



The value of P r \k =£ =£ k + 0.2 for k = to 4.0 is given in Table 1 for both the normal p.d.f. 

and f{x). 



We will now operate as if 



A.)-^-(-vi(^)r 



*We are indebted to Mr. Alan Kaplan of the AMC Inventory Research Office who suggested that we reconsider our values 
for a and b. 

**The absolute value of — ■ — is omitted. 



EOQ 



245 



Table 1 
Comparison o/f(x) to the Normal p.d.f. 



K 


A" 


B" 


C n (xY 


C f (xY 





0.0793 


0.1232 


0.0793 


0.1232 


0.2 


0.0761 


0.0928 


0.1554 


0.2160 


0.4 


0.0703 


0.0700 


0.2257 


0.2860 


0.6 


0.0624 


0.0527 


0.2881 


0.3387 


0.8 


0.0532 


0.0397 


0.3413 


0.3784 


1.0 


0.0436 


0.0299 


0.3849 


0.4084 


1.2 


0.0343 


0.0226 


0.4192 


0.4310 


1.4 


0.0260 


0.0170 


0.4452 


0.4480 


1.6 


0.0189 


0.0128 


0.4641 


0.4610 


1.8 


0.0132 


0.0097 


0.4773 


0.4704 


2.0 


0.0088 


0.0073 


0.4861 


0.4777 


2.2 


0.0057 


0.0055 


0.4918 


0.4832 


2.4 


0.0035 


0.0041 


0.4953 


0.4874 


2.6 


0.0021 


0.0031 


0.4974 


0.4905 


2.8 


0.0013 


0.0024 


0.4987 


0.4928 


3.0 


0.0006 


0.0018 


0.4993 


0.4946 


3.2 


0.0004 


0.0013 


0.4997 


0.4959 


3.4 


0.0001 


0.0010 


0.4998 


0.4969 


3.6 


0.0001 


0.0008 


0.4999 


0.4977 


3.8 


0.0001 


0.0006 


0.5000 


0.4982 


4.0 


0.0000 


0.0004 


0.5000 


0.4987 



— -=S/t + 0.2 \fom(x)=— i=exp 



•V2^ 



\ 2<r 2 / 



■B = Pr\k* ^k+0.2\forf(x)=— c .exp(-V2 

cr J 2cr \ cr 



x — fj. 



x-V- 



} 



X-/A 



C(x)=Pr =£ s Ar-H0.2l for n(x) and /(*), respectively. 



246 V. J. PRESUTTI, Jr. AND R. C. TREPP 

Notice that this is strongly biased against negative k. In fact, by operating as if the probability dis- 
tribution could be described by Eq. (4), we are, in effect, saying that limit Pr[x > p. + kcr] = °°, a 

k— *— °° 

somewhat preposterous allegation. Nonetheless, this bias is considered desirable since in actual 
practice, we restrict k 3 s 0.* 

Expected Backorders 

In order to determine the expected number of units backordered per year (Bx) and the average 
number of units in a backorder status at a random point in time (Bt) , we need g(y) —the p.d.f. for 
the number of backorders at a random point in time. The following technique for obtaining g(y) was 
developed by Hadley and Whitin [1]. 

g(y)dy— probability that the number of units backordered is between y and y + dy, and 



£-(-Vl(«±5=£))* 



= the probability that the number of units backordered at time, t, is between y, y+ dy given that 
the inventory position (on-hand plus due-in minus due-out) is between p, p + dp at time t — t. But the 
probability that the inventory position is between p, p + dp is dpIQ. 

Therefore, 

^r^e, P (-V2(^)), r * 
and 

and finally 

(7) g(y) = ^ exp (- V2 ( Z ^)) d " exp( -V2QI<r)). 

Consequently, the probability of a stockout (P ou t) is 

(8) | o " g(y)dy = (^|) | (1 - exp (- V2QI*)) exp (- V2k). 

Therefore, 

(9) B x = DP 0Ut = (~^j -Q- (1 - exp (- V2QI<t) ) exp (- V2k) 

(seeRef. [1], page 178) 
and 

(10) B T = f~ yg{y)dy = -^-^j (1-exp (- V2QI<r)) exp (- V2A). 

Other Inputs 

We must now obtain the average number of procurements per year, the expected inventory 
position (E[IP]), and the expected on-hand assets (E[OH]). Clearly, the average number of procure- 
ments is D/Q. 



*The nonnegativity requirement for A: is a result of the desire to make this model acceptable to practicing logisticians. 



EOQ 



247 



It is also well known that 



and 

By definition, 

Therefore, 



£[Due-in] = \x (see Ref. [1], page 187) 
E[IP]=r+QI2 = ix + ka + QI2. 

E[IP] =E[OH] +£[Due-in] -E[B/0]. 

E[OH] =E[IP] -£[Due-in] +E[B/0] 
= fi + ko- + QI2-ix + B T 
= k<r + QI2 + B T . 

FOUR INVENTORY MODELS 

We are now ready to proceed to the four models which will enable us to determine Q and r for each 
item in such a way as to minimize the holding cost (applied either to E[OH] or E[IP]) plus the pro- 
curement cost, subject to an overall backorder constraint. This constraint may be either the total units 
backordered per year (i.e., the stockout penalty is independent of the duration of the stockout) or the 
average number of units in a backorder position (i.e., the stockout penalty depends on the duration of 
the stockout). In either case, the backorders may also be weighted by item essentiality. 

MODEL I: Stockout penalty independent of duration of stockout; holding cost applied to on-hand 
inventory. We wish to minimize 



(LI) 

subject to 

(1.2) 



A AiD t - 

Z-nT + Z aiCi 

i=l v ' t=l 



kcri+y+^^(l-exp(-V2^) 



A 0.5 ZiDjO-j 



(l-exp(-V2^) exp (- V2ki) *S t, 



£V2 Q t 

where Z, represents item essentiality. The Method of Lagrange yields 

QidiCi 



(1.3) 

(1.4) 



*' = "V? ln 



-V^))(^p-XZ (J D, 



2 (afim -=V2-7,,and 



VJ 



-kZiDi 



(1.5) o = -^ + ^i 
Qf 2 



aidCTi 



V2" 1-exp [-V2 



£)) 



1-exp [-V2 



a-iJJ V2 



Qi 



<Ti 



exp 



(-^5) 



For reasons which will soon become evident, we will not solve (1.5) for Q t until later. 

MODEL II: Stockout penalty independent of duration of stockout; holding cost applied to in- 
ventory position. 



"AM" ( Q t 

Z ~q~ + Zj aiCi[fJ<i + ki(Ti + Y 



248 V. J. PRESUTTI, Jr. AND R. C. TREPP 

We wish to minimize 

(III) 

subject to 

" 0.5 ZiDm I ( _ Qi\ , 

The Method of Lagrange yields 

1 



(113) 



"v^ ln 



QidjCj 

0.5 (-A) Z,Z) ( fl-exp (- V2" 



<?.■ 



(H.4) 



A = — — — V cnatCi, and 
V2 Ti £i 



i" 



(115) = -^i + ^- 



aiCicn 



Vf(l-exp(-V2a)) 



1 — exp 



-V2& 






MODEL III: Stockout penalty depends on duration of stockout; holding cost applied to on-hand 
inventory. 

We wish to minimize 



(IH.l) 
subject to 
(III.2) 



|)^+|;aiC i [fe(ri+^+^g(l-exp(-V2^))exp(-V2A i )] 



i^^(l-ex P (-V2^))ex P (-V2fe)^ i 8. 



The Method of Lagrange yields 



(III.3) 



kt — ;= In 

V2 



V2(W, 



0.5(7, ( 1 -exp (- V2 & ) ) (aid-kZi) 



(III.4) 



1,(5^=^. and 



(III.5) °=-^ +S f- 



aiCicn 



V2(l 



exp 



V2 



CTi 



1 — exp 



V2& 



<?, 



-^expf-V^ 



o-j 



CXi 



MODEL IV: Stockout penalty depends on duration of stockout; holding cost applied to inventory 
position. 



EOQ 



249 



We wish to minimize 
(IV.l) 

subject to 

(IV.2) 



y -jt + i aiCi ( /*■• + kiai + 9 

i=l v* i = 1 \ ^ 



The Method of Lagrange yields 



(IV.3) 



ki= ~vz ln 



V2 QidiCi 



0.5 (- \)ZiO-, 1 -exp - V2 



(J, 



(IV.4) 



(IV.5) 



" cr ( aiCi 
: > — ;^t:, and 



<?? 2 " 



atCjCTi 



V? (l-exp - V2" 



-a 



<Ti 



l-exp l-v^ 



5i 



^exp(- 



(3i 



cr, 



Notice that in all four models Qi is the same (Eqs. (1.5), (II. 5), (HI. 5), and (IV.5)). The quantity 
Qi is completely independent of the backorder constraint regardless of whether or not the constraint 
is time-weighted. We also note that we cannot solve explicitly for Q. Conceptually, there is no problem. 
The exact solution for Q can be obtained through the use of an iterative routine (assuming access 
to a computer). An alternate and perhaps a more suitable approach is to recognize that in most cases 



(!-- p(-Vjfi)) 



] 



and can therefore be ignored. Consequently, 

AiDi 



(13) 
and 

(14) 



atCi ajCiCTi 
~ ~ Qf + 2 V2C, 



Qt V2 + V a iCi + lv2T 



Notice that as <t, — > 0, Eq. (14) approaches the well known Wilson lot size formula. 

DISCUSSION -MODEL APPLICATION 

Holding Cost Considerations 

We wish now to discuss the difference between applying holding cost to the on-hand inventory 
(Models I and III) and applying it to the inventory position (Models II and IV). In spite of tradition, we 
feel it is more logical to apply the holding cost to the inventory position (see Ref. [4] for an argument 
supporting this point of view). From a practical point of view, the difference in results is small; however, 
the inventory position models have two very desirable attributes. First, there is an explicit solution 



*It is easily shown that this value is an upper bound for the exact solution. 



250 



V. J. PRESUTTI, Jr. AND R. C. TREPP 



for X, the Lagrangian, in these models — which is not the case when the holding cost is applied to the 
on-hand inventory. Second, in the inventory position models the safety level expressions are equivalent 
to those that would have been obtained had we somehow arrived at the buy quantities ((Vs) and then 
independently minimized the cost of holding safety stock subject to the same backorder constraints 
used in Models I through IV. That is, under the inventory position models, the inventory manager 
who minimizes backorders subject to a safety level cost constraint would arrive at A:,'s equal to the A;,'s 
that he would have calculated had he minimized the total costs of ordering and holding inventory subject 
to a backorder constraint. This is especially important when the optimal (Vs are not used: for example, 
when the policy is to purchase 3, 6, 9, or 12 months' worth of stock, or even when Qi is determined by 
the Wilson lot size formula. If the same holding cost factor is applied to all items (i.e., at = a for all t), 
the results obtained in models II and IV are equivalent to those determined when minimizing the safety 
level investment subject to a backorder constraint. 

Negative k's 

Although all four models are strongly biased against negative &'s, as long as there is a heterogeneous 
inventory with some high cost, low demand items, there will be some items for which the models will 
specify negative &'s. It could be argued that if even with a strong bias negative A:'s are obtained they 
should be accepted. We have no philosophical quarrel with this argument; however, in this paper, we 
insist that there be no negative safety levels. In order to accomplish this, we will set all negative A's 
equal to zero. Consequently, we will have more protection than the models call for and thus have fewer 
backorders, but a higher cost. The problem is to determine the X that will allow, taking into account the 
nonnegativity of k, an acceptable number of backorders at a reasonable cost. Since we believe that 
Model IV is most appropriate for the Air Force's nonreparable items, we will use it to describe the 
methodology for determining A. 

Employing a simulation using actual data for about 5,000 Air Force items (for which seven years 
of history is available) we will select a range of values for [—A.]. Each value of [—A] will yield 5,000 
A;'s — one for each item. The nonnegative A:'s will be accepted and the negative &'s will be set equal to 
zero. We will then, using these A;'s which are now all nonnegative, obtain the average number of 
essentiality-weighted units in a backorder status divided by the number of essentiality-weighted units 
demanded annually (call this ratio /?«), and the cost of safety stock divided by the dollar value of annual 
demands (call this ratio R t ). Subsequently, Figures 1 and 2 will be plotted as shown below: 



-3t 




Given Figures 1 and 2, the inventory managers can select either R B or R$ (hopefully they will make 
the selection by considering both the cost and the effectiveness aspects) and then choose from the 
appropriate graph the value of [— X]. The ability to choose on the basis of fund availability is extremely 
important since, in practice, the inventory manager is considerably more likely to be subject to a 
system constraint of dollars available for safety level investment (out of pocket cost) than to a back- 



eoq 251 

order constraint. Hadley and Whitin emphasize this point in their* discussion of budget constraints. 
They conclude: 

"None of the inventory models discussed in operations research literature appears to provide 
an adequate treatment of budgetary constraint. In the first place, budgets are usually set in 
terms of actual dollars to be spent, whereas the costs included in inventory control models 
are often not the out-of-pocket cost type."* 
An example of how the inventory manager (IM) might use these graphs (possibly to reduce "actual 
dollars to be spent") is as follows: 

Suppose the IM submits a budget request, for operating and safety stock, which allows 20 percent 
of the dollar value of annual demands for safety stock. If he receives the monies requested he can 
obtain from Figure 2, the [— A.] that will yield an Rt of 20 percent. If he gets M dollars less than 
requested, and that M dollars represents 5 percent of annual dollar demands, he can set R% 
equal to 15 percent (20% — 5%) and obtain the new value for [— \]. 
While this paper doesn't concern itself (other than the portion immediately preceding) with 
constrained optimization, we have applied the demand distribution described in this paper to models 
that do constrain the buy quantities (i.e., the (?'s), the on-hand inventory, and the inventory position 
as well as the safety stock. None of these methods, however, is satisfactory when it comes to reducing 
actual dollars to be spent by very large amounts below computed requirements. All of these constraining 
models reduce dollars spent by reducing the reorder level (i.e.. by reducing k) and/or by reducing the 
order quantities in some manner. We cannot, however, reduce the reorder level greatly without paying 
a very large price in terms of backorders. If most of our expenditures were for items with multi-year 
buy quantities, then we could temporarily reduce actual dollars to be spent by reducing these buys 
to a year or less. Unfortunately, by far the greatest majority of dollars are spent on items for which the 
buy quantity is on the order of three months worth of stock — which theoretically eliminates the pos- 
sibility of obtaining a large reduction in actual dollars to be spent by reducing the buy quantity. 

Backorder Costs Known 

The reader will note that in the event the stockout cost (s,) is known, the problem is simplified. 
The (?,'s remain unchanged and the A,'s of each of the models have a 5; in place of [— A] Z, (see Eqs. 
(1.3), (II.3), (III.3), and (IV.3)). There is no [- A] . 

Final Note 

The Air Force is planning to test all four models, nevertheless, at the present time, it favors (un- 
officially) Model IV. This model has also been brought to the attention of the Department of Defense 
by the Air Force for consideration during a series of recent joint services meetings designed to further 
standardize inventory policies and objectives among all DoD components. 

REFERENCES 

[1] Hadley, G. and T. M. Whitin, Analysis of Inventory Systems (Prentice Hall, Englewood Cliffs, New Jersey, 1963). 

[2] Hadley, G. and T. M. Whitin, "A Review of Alternative Approaches to Inventory Theory," The RAND Corporation, Santa 

Monica. California. RM-4185-PR (Sept. 1964). 
(3| Parker, L. L.. "Economical Reorder Quantities and Reorder Points with Uncertain Demand," Nav. Res. Log. Quart.. //, 

351-358 (Dec. 1964). 
[4] Schrady, D. A., "Continuous Review Inventory Practice in the Navy Application," Report #NPS55o9031A, United States 

Naval Postgraduate School, Monterey, California (Mar. 1969). 



* Reference [2|, page 152. 



THE RELATIONSHIP BETWEEN DECISION VARIABLES AND PENALTY 
COST PARAMETER IN (Q, R) INVENTORY MODELS 

Alan J. Kaplan 

AMC Inventory Research Office 

Army Logistics Management Center 

Fort Lee, Virginia 



ABSTRACT 

This paper is concerned with the optimum decision variables found using order quantity, 
reorder point (Q, R) inventory models. It examines whether the optimum variables 
(Q* and R*) are necessarily monotonic functions of the backorder cost parameter (or 
equivalently of the performance objective). For a general class of models it is proved that 
R* must increase as the performance objective is raised, and an inequality condition is 
derived which governs how Q* will change. Probability distributions of lead time demand 
are cited or found for which Q* increases, Q* decreases, and Q* is independent of increases 
in performance objectives or backorder cost parameter. 

INTRODUCTION 

This paper is concerned with the optimum decision variables found by using order quantity, 
reorder point (Q, R) inventory models. It examines whether they are necessarily monotonic functions 
of the penalty (backorder) cost parameter (or equivalently of the performance objective). It attempts 
to make general statements about the direction of such monotonicity. Such information is of particular 
interest in contexts where the true cost of a backorder is difficult to estimate — and the performance 
objective is determined in large part by its impact on inventory levels. 

It is easy to be misled by intuition. Our own empirical experience was with a (Q,R) model in which 
backorder cost is proportional to the length of time a backorder is outstanding, and lead time demand 
is normally distributed. For this case, the optimum reorder point (R*) always increased as the penalty 
cost parameter (P) increased, while the optimum order quantity (Q*) decreased. This is illustrated 
in Herron's [2] dimensionless graphs for the normal distribution. The conclusion can be proved analyti- 
cally for deterministic demand and penalty proportional to time on backorder by solving directly for Q *; 
e.g., see Hadley and Whitin's [1] solution. 

On the other hand, Victor Presutti [3], found Q* to be independent of the penalty cost parameter 
when he used a pseudo-exponential distribution for lead time demand. t That finding prompted this 
paper. 

We initially assume penalty cost is proportional to the length of time a backorder is outstanding 
and begin by reviewing the relevant aspects of the appropriate (Q, R) model. Next we obtain some 
useful theoretical conclusions, including a proof that R * must always increase when P does and an 
inequality condition which governs how Q* will behave. We apply these results, and then in the last 
section extend them to other {Q, R) models. 



t/(x) = A.ie x " forx^O, and/(jc) denned elsewhere so that / f(x)dx=\. Presutti's work is based on Lawrence Parker's 
approximation to the normal complementary cumulative. 

253 



254 A - J- KAPLAN 

The theoretical conclusions require the use of an approximation for average backorders (Eq. (2)). 
This approximation is commonly used in the actual computation of /?* and Q* [1], [2J, and does not 
negate the value of the results. 

THE MODEL 

The model is the basic steady state, continuous review (Q, R) model applicable when penalty 
cost is proportional to the length of time a backorder. is outstanding. For the model to give strictly 
optimal answers in the general case, procurement lead times must be constant and demands must 
be for one unit at a time. However, the model is being used profitably even when these conditions do 
not hold. For ease of exposition we will assume that the lead time demand distribution is continuous. 
One good reference, particularly for providing the basis for Eqs. (1) and (2), is Hadley and Whitin [1]. 

Let 

TC = average total cost per year; 
H= holding cost per year per unit of stock on hand; 
P= penalty cost per year per unit of stock on backorder; 
A = administrative and/or setup cost; 
A — average yearly demand; 

A/. = expected demand during procurement lead time; 
B — average amount on backorder at a random point in time; 
G(x) — probability that lead time demand equals or exceeds x; 
L = lead time. 
Then 

TC(R,Q)=A(klQ)+H- (R + QI2-\,. + B) +P • B = A • {KIQ)+H- (R + QI2-k,.)+ P'B 

(1) 

P' = P + H 

B = ±r (x-R)[G(x)-G(x + Q)]dx. 
v Jk 

For computational purposes, the following approximation is often used for B: 

(2 ) B = jjf~(x-R)G(x)dx. 

The approximation is valid to the extent that the probability of demand exceeding R + Q is small. 
We will utilize this expression for B. 

THEORETICAL RESULTS 

Let 

Qi, R, be the optimal decision variables for P'=/ > t ; 
Qi, Ri be the optimal decision variables for P' =P 2 ; 

Pz>Pr, 
Bt be expected average backorders for decision variables Qi, Rr, 



F(Q,R) 



(Q,R) INVENTORY MODELS 255 

dB/dQ 
SB/dR' 



The following two theorems will be proved: 

THEOREM 1: R* (the optimal R) must increase as P increases. 

THEOREM 2: The sign ofdF(Q, R)/dR, if it is consistent, determines whether ()* will also increase. 

The key to these results is the corollary to Lemma 1 that the algebraic sign of Q 2 — Q\ is determined 
by the algebraic sign of F(Q 2 , R 2 ) — F(Q U R t ). Lemma 2 and its corollary demonstrate the intuitively 
"obvious": that expected backorders (B) must decrease in the optimum as P increases and, therefore, 
that it is not possible for both Q* and R* to decrease. Lemma 3 shows that dF{Q, R)ldQ is always 
negative. Lemma 2 is used in proving Theorem 1, while Lemma 3 is used in proving Theorem 2. 

LEMMA 1: The algebraic sign of F(Q 2 , Rz) — F{Q\, R\) is the same as the algebraic sign of 
Qi~Q\- (A zero value for one implies a zero value for the other.) 

PROOF: Taking the partial derivatives of the total cost equation, Eq. (1), with respect to Q* and 
/?*, and setting the derivatives to zero gives two conditions: 

(3a) (dBldQ*)-P'-^-AI(Q*y + H/2 = 0; 

(3b) (dB/dR*)-P' + H=0. 

By algebraic manipulation (transposing and dividing one equation by the other), this gives 

F(Q*,R*) = [A -AI(Q*) 2 -HI2] +-H. 
Hence 

F(Q 2 , R 2 ) -F(Q U /?,) = [X-AKQ.r-X-AKQ.y] --//. 

If Q 2 exceeds Qi(i.e., Q2 — Q1 is positive), then both the numerator and denominator on the right 
hand side of Eq. (4) are negative; F(Q 2 , R 2 ) — F(Qi, Ri) is then positive. By similar reasoning the 
rest of Lemma 1 is shown to hold. 

COROLLARY: The algebraic sign of Qt — Q\ is the same as the algebraic sign of F(Q 2 , R 2 ) 
— F(Q\, Ri). A zero value for one implies a zero value for the other. 

LEMMA 2: The expected value of backorders (in the optimum) must decrease as P increases. 

PROOF: To simplify the explanation let Oi stand for operating costs corresponding to decision 
variables (Q u Ri) so that Total Cost = O t + B t ■ P' . 

We use proof by contradiction. We initially assume B 2 > B t . By the optimality of Q u Ri for/ > '=P 1 , 

Using P, > P, and B>>Bi, we find 

O. + BrP^BriPi-Pt) > O. + fi.P. + fi.^Pi-P,). 
Hence 

2 + B 2 P 2 >0 1 + Bi P 2 . 

But this last inequality contradicts the optimality of Q 2 , R 2 when P' =P 2 . 

COROLLARY: Both Q * and R * cannot decrease as P increases (proved by Lemma 2 and inspection 
of Eq. (2)). 

LEMMA 3: dF{Q, R)/dQ is always negative. 

PROOF: Using Eq. (2), we get by calculus, 



256 A. J. KAPLAN 

(5a) % = ~7p\l ('-WW** 

dB 1 f- v , 
= ~Q in G(x)dX - 



and 
(5b) 

Hence, 



dR 



(6) F(Q, R) 



~j^(x-R)G(x)dx ±j^ (x-R)G(x)dx g 

— I™ G{x)dx !" G(x)dx [ X G{x)dx 

(JjR JR JR 



The denominator does not depend on Q. Clearly, the numerator decreases as Q increases. 
THEOREM 1: R* increases as P increases. 
PROOF: Using Eq. 6, we have 

(7a) F(Q 2 ,R 2 )=B 2 +\G(x)dx 

J «2 

and 

(7b) F(Q u R 1 )=Bi^r('G(x)dx. 

JRi 

Suppose R 2 < Ri. Then the denominator of (7a) will be larger than the denominator of (7b). By 
Lemma 2, B 2 <B\, so this means F{Q 2 , R 2 ) <F(Qi, Ri). By the corollary to Lemma 1 this implies 
Qz<Q\. By the corollary to Lemma 2, however, it is not possible for Q 2 < Q\ if R 2 < Ri since P 2 > Pi. 
Proof by contradiction. 

THEOREM 2A: If dF(Q, R)/dR is consistently negative, then Q* must decrease as P increases. 

Since R 2 >Ru and dF(Q, R)/dR is given as negative, F(Q U R 2 ) <F(Q l , R^. Now if Q 2 ><?,, 
then by Lemma 3 this means F(Q 2 ,R 2 ) < F(Qi, R 2 ). By combining the two inequalities, we get 
F(Q 2 , R 2 ) < F(Qi, Ri). However, by the Corollary to Lemma 1, this implies Q\>Q 2 . Proof by 
contradiction. 

THEOREM 2B: If dF(Q, R)/dR is consistently positive, then Q* increases with P. 

The proof is analagous to that of Theorem 2 A. Since R 2 > Ri, and dF(Q, R)ldR is given as positive, 
F(Q U R 2 ) >F(Q U /?,). Now if Q 2 <Q U then by Lemma 3 this means F{Q 2 ,R 2 ) >F(Q U R 2 ). By 
combining the two inequalities we get F(Q 2 , R 2 ) > F(Q t , R\). But by the Corollary to Lemma 1 this 
implies @i < (?2. Proof by contradiction. 

THEOREM 2C: If F{Q, R) is independent of R, then Q* is independent of P. For this case 
F(Q U R 2 )=F(Qi, Rt). If <? 2 < <?i then by Lemma 3, F((? 2 , /? 2 ) > F{Q X ,/?,), which by the Corollary 
to Lemma 1 implies Q 2 >Q\. Similar reasoning shows that the assumption ()i < Q 2 also leads to a 
contradiction. 

Applications 

There is a strong temptation to reason intuitively that Q* must increase as penalty cost increases: 
higher reorder quantities reduce the probability of stock out; this advantage becomes more significant 
as the penalty cost parameter is increased, while the impact of Q on ordering and holding costs is 
unchanged. Why is this reasoning wrong? 

The flaw is that when penalty cost increases, it becomes worthwhile to increase /?*, so that chang- 
ing Q* has less impact on backorders than before, and therefore may have less impact on penalty 



(Q,R) INVENTORY MODELS 257 

costs (the product of backorders and the penalty cost parameter). The decisive factor relates to the 
relative effectiveness in reducing backorders of augmenting Q* vs augmenting R*. As R increases, 
does the relative effectiveness of augmenting Q increase? Only then will Q* increase as P increases. 

It is difficult to apply Theorem 2 analytically to most probability distributions, i.e., to determine 
analytically the sign (+, — , or mixed) of dF(Q, R)dR. No attempt was made to carry on an extensive 
computer aided mapping of sign vs probability distribution 

We did attempt to find a case when Q* would increase as penalty cost increases. Recall that for 
the normal distribution Q* decreases and that for exponential type distributions Q* is unchanged. 

Let the Weibull distribution be defined: 

f(x) = aXx"- 1 exp (- Xx a ); x^0,Xa^0. 

Then for a= V2, we have analytically shown that dF(Q, R)/dR is positive for R 3= 0, i.e., this is a 
case where Q* increases with R. Empirical work indicates that this conclusion about the Weibull 
distribution also holds true when the precise expression for average backorders (Eq. (1)) is used.t 

For a=l the Weibull reduces to the exponential. It is simple to show dF(Q, R)/dR — for the 
exponential regardless of which expression for B is used. 

As a matter of interest, the Weibull with a= V2 was originally examined because of an inference 
made from an analogy between determination of the sign of dF(Q, R)/dR, and the sign of the deriva- 
tive of the failure rate function of reliability theory, the latter being a standard problem. 

EXTENSIONS 

The general linear form for penalty cost is P + P • t where t is the time on backorder. This trans- 
lates into total expected penalty cost =P ■ B + P ■ B or, in words (cost per unit backordered per year 
on backorder) X (average amount on backorder) + (fixed cost per backorder occasion) X (number of 
demands backordered each year). That part of the total cost equation, Eq. (1), dependent on B or B 
becomes P' B + P ■ B where P' = P + H. 

The theory derived in this paper still holds provided that as the penalty cost parameters are 
increased, P remains equal to k • P' for constant k (or P — 0, the special case already considered). 
If the cost item H • B is left out of the total cost equation, Eq. (1), the condition is simply that P — k • P, 
which is more natural. H • B is normally very small, and in fact does not appear as a cost if holding 
cost depends on inventory position (including net due-in) rather than on hand. 

An approximate expression must be used for B as well as B. 1 1 

Let 



fjl'W 



B T =kB + B. 

For F T {Q, R), defined in terms of dBrldQ and BB T lBR, we have 

Br 



Ft(Q,R) = 



(*G(x)dx + k-k-G(R) 



tMiss Chung-Mei Ho did much of the analytical investigation on the application of Theorem 2, and Mr. Robert Deemer 
did the empirical work for the Weibull. 

ttHadley and Whitin, chapter 4 has a derivation. 



258 A. J. KAPLAN 

By inspection, all three lemmas proved hold with the substitutions of Bt for B, P' • Bt for P'B, and 
F T (Q, R) for F(Q, R). Therefore, the two theorems hold with these substitutions. 
Using the precise expression for B causes difficulties even when P is 0. Then 

B-T [y-(R + Q)]g(y)dy , C( . 

F{Q,R)=—±** , where g ( y ) =-&&-. 

[G(x)-G(x + Q)]dx . y 

It is no longer apparent that Lemma 3 must always hold, or that in the proof of Theorem 1 that 
the numeraton F(Q 2 , Rz) < numerator a£F{Qu Ri) if Qi > Q\- 

REFERENCES 

[1] Hadley, G. and T. M. Whitin, Analysis of Inventory Systems (Prentice Hall, Englewood Cliffs, N.J., 1963). 

[2] Herron, David P., "Use of Dimensionless Ratios to Determine Minimum-Cost Inventory Quantities," Nav. Res. Log. Quart. 

13, 167-176(1966). 
[3] Presutti, V. J. and R. C. Trepp "More Ado about Economic Order Quantities (EOQ)", Nav. Res. Log. Quart. 17, 242-250 

(1970). 



NEWS AND MEMORANDA 

PROFESSIONAL DESIGNATION IN LOGISTICS MANAGEMENT 

The School of Systems and Logistics, Air Force Institute of Technology, in cooperation with 
the Society of Logistics Engineers (SOLE), has announced the beginning of a program of professional 
designation for qualified military and Civil Service logistics managers of the Department of Defense. 
Those who meet the criteria may be awarded an appropriate certificate of "Professional Designation 
in Logistics Management." Certification is not restricted only to members of the Society. 

The purpose of the program is to recognize the efforts of the individual to increase his educa- 
tional background in the broad world of logistics. It recognizes that true logistics management capa- 
bility must encompass more than expertise in one of the several disciplines which comprise the logistics 
spectrum. The idea is that logistics is inherently an integrating process of these several disciplines 
and that the effective logistics manager will be one who has knowledge in, and understands the view- 
point of, the various sub-functions which the separate disciplines represent. The professional stature 
of the functional specialists is recognized, but this program is aimed at the need lor a portion of these 
managers to become logistics generalists to perform the essential function of integrating the efforts 
of many functional specialists. 

Qualification for certification rests upon the types and forms of education and experience which 
contribute to the overall concept of integrative management in logistics. The courses of study outlined 
in the qualification listing were tailored with this concept in mind. The objective is to assure that those 
who qualify have progressed from one of the functional specialties to a broader realization of logistics 
management needs through exposure to education in related and supporting disciplines and functional 
areas. Thus, those who qualify for certification will have demonstrated a high degree of initiative in 
acquiring a broad educational and experience background. They will be considered to have demon- 
strated their competency for advancement to higher levels of responsibility in logistics management 
where their broadened point of view would aid decision-making of more cost-effective nature. 

The certification program provides for two levels of certificate award. The first is "Professional 
Designation in Logistics Management" and the second is "Professional Designation in Logistics 
Management — Advanced." The qualification criteria includes courses in five categories as indicated: 

\ Certificate Advanced 

Category A Logistics Integration 2 2 

Category Bl Materiel Acquisition 1 2 

Category B2 Materiel Distribution 2 2 

Category B3 Materiel Maintenance 1 2 

Category C Management Techniques 2 2 



Total Required 8 10 

Some allowance is made in the program for waiving certain educational course requirements on 
the basis of field experience. Such waiver, for a maximum of two courses, must be substantiated by 
the applicant. The determination of whether or not to accept the experience equivalency is left to the 
discretion of the certification committee of the School of Systems and Logistics. Applicants are urged 

259 



260 NEWS AND MEMORANDA 

to explain their request for equivalency and waiver by an adequate explanation of the experience 
(where, for how long, in what job, at what level, and so on). 

It is also likely that many will have taken higher level educational courses than those listed in 
the certificate categories. These may be accepted in lieu of the described courses if the applicant 
includes proof of completion along with a brief outline of material covered. 

Certificates will be approved by the School of Systems and Logistics and the SOLE National 
Education Committee through a certification board process. Applicants, who need not be members 
of SOLE but who must be members of DOD, are encouraged to submit their requests for certification 
to the School of Systems and Logistics (AFIT-SL), W-PAFB, Ohio 45433. Each request must be 
accompanied by the necessary back-up data. 

CHANGE IN OFFICE OF NAVAL RESEARCH ADDRESS 

Effective with the relocation of Department of the Navy activities from the Main Navy Building 
in Washington, the address for the Office of Naval Research, including the Naval Research Logistics 
Quarterly, is Arlington, Virginia. 

Correspondence for the Naval Research Logistics Quarterly should be addressed as follows: 
Managing Editor 

Naval Research Logistics Quarterly 
Office of Naval Research 
Arlington, Virginia 22217 



■fr U.S. GOVERNMENT PRINTING OFFICE i 1970 OL-386-049 



INFORMATION FOR CONTRIBUTORS 

The NAVAL RESEARCH LOGISTICS QUARTERLY is devoted to the dissemination of 
scientific information in logistics and will publish research and expository papers, including those 
in certain areas of mathematics, statistics, and economics, relevant to the over-all effort to improve 
the efficiency and effectiveness of logistics operations. 

Manuscripts and other items for publication should be sent to The Managing Editor, NAVAL 
RESEARCH LOGISTICS QUARTERLY, Office of Naval Research, Arlington, Va. 22217. 
Each manuscript which is considered to be suitable material tor the QUARTERLY is sent to one 
or more referees. 

Manuscripts submitted for publication should be typewritten, double-spaced, and the author 
should retain a copy. Refereeing may be expedited if an extra copy of the manuscript is submitted 
with the original. 

A short abstract (not over 400 words) should accompany each manuscript. This will appear 
at the head of the published paper in the QUARTERLY. 

There is no authorization for compensation to authors for papers which have been accepted 
for publication. Authors will receive 250 reprints of their published papers. 

Readers are invited to submit to the Managing Editor items of general interest in the field 
of logistics, for possible publication in the NEWS AND MEMORANDA or NOTES sections 
of the QUARTERLY. 



NAVAL RESEARCH JUNE 1970 

LOGISTICS VOL. 17, NO. 2 

QUARTERLY NAVSO P-1278 



CONTENTS 

ARTICLES 

Page 
Markov Chain Analysis of a Situation Where Cannibalization is the Only Repair Activity by A. J. 
Rolre 

Some Estimates of Reliability Using Interference Theory by M. Mazumdar 

The Optimal Burn-In Testing of Repairable Equipment by J. M. Cozzolino 

Goal Programming in Econometrics by W. A. Spivey and H. Tamura 

i 

Production and Employment Scheduling in Multistage Production Systems by C. Haehling von 
Lanzenauer 

A Mutual Primal-Dual Linear Programming Algorithm by M. Y. Harris 

Interception in a Network by R. D. Wollmer 

The Fixed Charge Problem by D. I. Steinberg 

On Computing the Expected Discounted Return in a Markov Chain by D. F. Hitchcock and J. B. 
MacQueen 

More Ado About Economic Order Quantities (EOQ) by V. J. Presutti, Jr. and R. C. Trepp 

The Relationship Between Decision Variables and Penalty Cost Parameter in (Q,R) Inventory 

Models by A J. Kaplan 

News and Memoranda 



OFFICE OF NAVAL RESEARCH 

Arlington, Va. 22217